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ABSTRACT 


Sr, Nd and Pb isotopic analyses of 477 samples representing 30 islands or 
island groups, 3 seamounts or seamount chains, 2 oceanic ridges and 1 oceanic 
plateau [for a total of 36 geographic features] are compiled to form a 
comprehensive oceanic island basalt [OIB] data set. These samples are 
supplemented by 90 selected mid-ocean ridge basalt [MORB] samples to give 
adequate representation to MORB as an oceanic basalt end-member. This 
comprehensive data set is used to infer information about the Earth's mantle. 
Principal component analysis of the OIB+MORB data set shows that the first 
three principal components account for 97.5% of the variance of the data. Thus, 
only four mantle end-member components [EMI, EMIT, HIMU and DMM are 
required to completely encompass the range of known isotopic values. Each 
sample is expressed in terms of percentages of the four mantle components, 
assuming linear mixing. There is significant correlation between location and 
isotopic signature within geographic features, but not between them, so 
discrimination analysis of the viability of separating the oceanic islands into those 
lying inside and outside Hart's (1984, 1988) DUPAL belt is performed on the 
feature level and yields positive results. 


A "continuous layer model" is applied to the mantle component percentage 
data to solve for the spherical harmonic coefficients using approximation 
methods. Only the degrees 0-5 coefficients can be solved for since there are only 
36 features. The EMI and HIMU percentage data sets must be filtered to avoid 
aliasing. Due to the nature of the data, the coefficients must be solved for using 
singular value decomposition [SVD], versus the least squares method. The F-test 
provides an objective way to estimate the number of singular values to retain 
when solving with SVD. With respect to the behavior of geophysics control data 





sets, only the degree 2 spherical harmonic coefficients for the mantle components 
can be estimated with a reasonable level of confidence with this method. 


Applying a "delta-function model" removes the problem of aliasing and 
simplifies the spherical harmonic coefficient solutions from integration on the 
globe to summation over the geographic features due to the properties of delta- 
functions. With respect to the behavior of geophysics control data sets, at least 
the degree 2 spherical harmonic coefficients for the mantle components can be 
estimated with confidence, if not the degrees 3 and 4 as well. Delta-function 
model solutions are, to some extent, controlled by the nonuniform feature 
distribution, while the continuous layer model solutions are not. 


The mantle component amplitude spectra, for both models, show power at 
all degrees, with no one degree dominating. The DUPAL components (EMI, 
EMII and HIMUJ, for both models, correlate well with the degree 2 geoid, 
indicating a deep origin for the components since the degrees 2-3 geoid is 
inferred to result from topography at the core-mantle boundary. The DUPAL 
and DMM components, for both models, correlate well [and negatively] at degree 
3 with the velocity anomalies of the Clayton-Comer seismic tomography model 
in the 2500-2900 km depth range [immediately above the core-mantle boundary]. 
The EMII component correlates well [and positively] at degree 5 with the 
velocity anomalies of the Clayton-Comer model in the 700-1200 km depth range, 
indicating a subduction related origin. Similar positive correlations for the geoid 
in the upper lower mantle indicate that subducted slabs extend beyond the 670 
km seismic discontinuity and support a whole-mantle convection model. 


Thesis Supervisor: Dr. Stanley R. Hart 
Title: Senior Scientist, Woods Hole Oceanographic Institution 
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CHAPTER 1 


INTRODUCTION 


PREVIOUS WORK 

That the Earth's mantle is heterogeneous is no longer a subject of 
controversy among geochemists, but the composition, the location and the 
geometry of these heterogeneities is very much in question. Direct sampling is 
not an option for studying the chemistry of most of the mantle, so products of 
indirect sampling, such as oceanic island basalts [OIB's] and mid-ocean ridge 
basalts [MORB], are invaluable for revealing the nature of the inaccessible 
mantle. Though the OJB's may be contaminated by interactions with the 
lithosphere or may sample large vertical sections of the mantle, they still retain 
the signature of their original source. 

Using various statistical methods and models, previous workers have 
defined what they believe to be the number of mantle component end-members 
required to represent the variation in the oceanic mantle data [OIB+MORB]. 
Early on, Zindler et al. (1982) used factor analysis to evaluate the oceanic data in 
five dimensions. Their analysis indicated that the oceanic data define a plane [the 
"mantle plane"], described by the mixing of three chemically independent 
components, two undifferentiated or slightly enriched mantle components and 
one MORB-type or depleted mantle component. 

Other workers have chosen five атаја or components to represent the 
data. Using a series of two-dimensional isotopic plots, White (1985) divided the 
oceanic data into five distinct basalt groups [MORB, St. Helena, Kerguelen, 
Society, and Hawaii]. He concedes that the five groups may be end-members 
which mix to form intermediate compositions, but he believes that each group 


either represents a distinct, internally homogeneous reservoir or that each group 
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is composed of a number of isotopically similar reservoirs. Likewise, Li ef al. 
(1991) proposed fives extremes, using non-linear mapping: Atlantic MORB 
[DMM], St. Helena [HIMU], Walvis [EMI], Samoa [EMII] and Ds [EMMI]. Non- 
linear mapping approximately preserves the geometric structure of the data by 
maintaining interpoint distances. Four of the five extremes of Li et al. (1991) 
are based solidly on samples trends from islands, but the Ds extreme is based 
only on that one sample. More data is needed to substantiate their fifth extreme. 

By far the majority of analyses indicate the existence of four end-member 
components for the oceanic mantle data. Using two-dimensional plots, Zindler 
and Hart (1986) defined the following four end-member components: depleted 
MORB mantle [DMM], high U/Pb mantle [HIMU], and two enriched mantle 
components [EMI and EMII], with possibly two other components prevalent 
mantle composition [PREMA ] and bulk silicate Earth [BSE]. Eigenvector 
analyses by Allegre et al. (1987) agree with the four component model of 
Zindler and Hart (1986). The four extremes of Allčgre et al. (1987) are 
[correspond to]: extreme MORB [DMM|]; St. Helena, Tubuai and Mangar islands 
[HIMU]; Kerguelen, Gough, Tristan da Cunha and Raratonga islands [EMI]; and 
Sào Miguel and Atui islands [EMII]. Hart (1988), using an augmented data set 
and two-dimensional plots, concluded that the four end- members proposed by 
Zindler and Hart (1986) are valid representations of the extremes of the oceanic 
data. He resolves White's (1985) groupings into his own four component system 
as Dt nı: = Hart]: MORB = DMM, Society = EMII, St. Helena - HIMU, 
Hawai = EMI, with the suggestion that White's fifth group, Kerguelen, is a 
mixture of EMI and EMII. In addition, Li et al. (1991) also noted a tetrahedral 
structure to the data, when using factor analysis with varimax rotation, with the 
following four extremes: Atlantic MORB [DMM], Mangaia [HIMU], Samoa 
[EMI] and Walvis[EMI). 





E 


One scenario for the genesis of the three unusual mantle components is put 
forth by Hart (1988). He proposes that HIMU, enriched in U, is probably 
generated by intra-mantle metasomatism, that EMI corresponds to a slightly 
modified bulk-earth compositon and that EMII can be explained by the recycling 
of sediments during subduction. The proposed formation mechanisms in no way 
limit the geometry of the mantle needed to generate the heterogeneities and, as 
such, a wide variety of models have been proposed. A whole mantle convection 
model might portray the enriched mantle components as blobs floating around in 
a depleted mantle matrix (Zindler and Hart, 1986) or perhaps as an accumulated 
layer of subducted oceanic crust and sediment at the core-mantle boundary that 
reaches the surface in mantle plumes (Hofmann and White, 1982). A two-layer 
convection model might rely on a depleted upper mantle feeding the mid-ocean 
ridges and an enriched lower mantle feeding oceanic islands via mantle plumes 
(Dupre and Allegre, 1983) or require a depleted upper mantle, a primitive lower 
mantle and an accumulated layer of subducted oceanic crust and sediment at the 
670 km discontinuity that supplies the enriched components via mantle plumes 
(White, 1985; Allégre and Turcotte, 1985). Anderson (1985) even proposes a 
three-layer convective model with the geochemical contrasts occurring only in 
the upper mantle with a depleted lower part that supplies the mid-ocean ridges 
and an enriched upper part from subduction of oceanic crust and sediment. 

A deep origin for the enriched components is indicated by Hart's (1984) 
large-scale isotopic anomaly, the DUPAL anomaly, characterized by the 
concentration of the enriched mantle components in a band from 2° S to 60° S. 
Qualitatively, countours of the anomaly criteria [A7/4, A8/4 and ASr (Hart, 
1984)} correspond to long-wavelength [and thus deep] geophysical quantities 
(Hart, 1988). Other researchers oppose this deep origin interpretation, citing the 


nonuniform distribution of hotspots as the reason for the pattern (White, 1985) 
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or arguing that the DUPAL compositions occur in scattered locations and do not 
cover a coherent geographic area (Allegre et al., 1987). 

The purpose of this thesis is three-fold: (1) to address once again the issue 
of the number of mantle end-member components needed to represent the 
oceanic mantle data, (2) to statistically test the viability of the DUPAL distinction 
as a means of characterizing the OIB data and (3) to try to pinpoint the source 
depth of the ennched mantle components by expanding their relative abundances 
in spherical harmonics and comparing their expansions to those of known 


geophysical quantities. 


DATA 

The majority of this study focuses on Sr, Nd and Pb isotopic analyses of 
volcanic rocks from oceanic islands, seamounts, ridges, and plateaus. All of 
these geographic features overlie oceanic crust, with the exception of Nunivak 
Island on the Alaskan Continental Shelf, and none of them is directly associated 
with seafloor spreading, with the exception of Iceland, which has a mixture of 
mid-ocean ridge and hotspot influences. Essentially, the data set is that compiled 
by Zindler et al. (1982) and later augmented by Hart (1988), with some 
additional recent analyses (Appendix). Samples in the data set are mainly basalt, 
with some gabbros and trachybasalts; trachytes and other silica-rich rocks 
relative to basalt [roughly SiO > 50%] are excluded. The majority of the 
samples are of Cenozoic age, with the exception of the Walvis Ridge, Rio Grande 
Rise and New England Seamounts samples, with ages up to 100 Ma. If a choice 
is given, analyses of leached samples are preferred over analyses of unleached 
samples. In addition, only single samples for which there are Sr, Nd and Pb 
analyses are included. For consistency, Sr data is adjusted to 0.70800 [E& A 


standard] or 0.71022 {NBS SRM 987 standard] and Nd data is adjusted to 
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0.512640 [BCR-1 standard] or 0.511862 [La Jolla standard] or 0.511296 [Spex 
standard]. 

In this data set, referred to as the OIB data set, there are 477 samples 
representing 30 islands or island groups, 3 seamounts or seamount chains, 2 
aseismic oceanic ridges and 1 oceanic plateau (Figure 1.1 and Table 1.1). The 
isotopic means and standard deviations for the OIB data are listed in Table 1.2. 

Since MORB is considered to be one of the mantle component end- 
members (Zindler et al., 1982; White, 1985; Zindler and Hart, 1986; Hart, 
1988), any attempt to choose end-members should include MORB data. For this 
reason, a second data set is created using the OIB data and a selection of 90 
MORB samples (Appendix), the OIB+MORB data set (Table 1.3). The criteria 
for choosing OIB samples applies to the MORB samples as well. Isotopic means 


and standard deviations for the OIB+MORB data are listed in Table 1.2. 


ORGANIZATION 

The main thust of this work is to characterize the OIB data and to search 
for possible correlations between the geochemical signatures of OIB's and 
geophysical quantities, such as the geoid and seismic tomography, that might help 
pinpoint the depth[s] of the OIB reservoir[s]. 

Chapter 2 explores the nature of the OIB isotope data. With the help of 
principal component analysis, the data is expressed in terms of percentages of 
four mantle component end-members. Spatial correlation testing reveals the 
relationship between geographic distance from island to island and feature to 
feature and the "isotopic distance" between samples. Discrimination analysis, 
both nearest-neighbor and graphical, is used to test the viability of separating the 


oceanic islands into two groups, inside and outside the DUPAL belt. 
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Chapter 3 applies a "continuous layer model" to the mantle component 
data, as an assumed geometry for the OIB reservoir, in order to solve for the 
spherical harmonic coefficients. The problem of aliasing 1s addressed with the 
relationship of variation in mantle components to distance between features. 
Approximation methods are used to solve for the coefficients. Geophysical data 
sets are constructed, using GEM-L2 geoid coefficients, to serve as controls 
against which to judge the success of the approximation methods. 

Chapter 4 applies a "delta-function model" to the mantle component data 
to provide a mathematically more robust solution for the spherical harmonic 
coefficients. The delta-function approximation removes the problem of aliasing, 
but generates a solution dependent upon feature location. The same geophysical 
data sets are used again to judge the success of the delta-function approximation. 

Chapter 5 compares the mantle component spherical harmonic solutions 
for the two models in terms of their amplitude spectra, how well they correlate 
with the geoid, how they are affected by the nonuniform feature distribution and 
how well they correlate with the Clayton-Comer seismic tomography model. 
The implications of these results and recommendations for further research are 


discussed. 
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Table 1.1. Geographic features represented in the OTB data set, vvith their 
components, number of samples [in braces] and references indicated. 


Feature 


Ascension [5] 


Amsterdam/St. Paul [11] 


Azores [6] 


Balleny [3] 


Cameroon Line [18] 


Cape Verde Islands [41] 


Christmas [13] 


Cocos [3] 


Comores Archipelago [14] 
Cook-Austral Islands [26] 


Crozet Islands 19) 
Fernando de Noronha 116) 
Galapagos İslands 111) 


Gough 12) 


Hawaiian Islands [73] 


Components 


Amsterdam [5] 
St. Paul [6] 


Faial [1] 
São Miguel [5] 


Bioko [5] 
Pagalu [1] 
Principe [3] 
São Tomé [9] 


Fogo [6] 

Maio [9] 

Sao Antao [10] 
Sao Tiago [13] 
Sao Vincente [3] 


Aitutaki [4] 
Atu 16) 
Mangia [5] 
Mauke [3] 
Raratonga 18) 


Hawaii [14] 
Kahoolawe [13] 
Kauai [2] 

Lanai [4] 


References! 


7,34,35 


38 
38 


.— 


IE 1:25 
li; 

p 
28,32 


28 
37 





Table 1.1. Continued. 


Feature 


Hawaiian Islands [73] 


Iceland [7] 


Juan Fernandez Islands [4] 


Kerguelen Plateau [41] 


Louisville Seamount Chain [4] 
Marion/Prince Edward [4] 
Marquesas Archipelago [11] 


Mascareignes [8] 


New England Seamounts [6] 


Nunivak [2] 
Pitcairn [19] 
Ponape [1] 

sala Y Gomez [1] 
Samoa İslands 1341 


San Felix/San Ambrosio [5] 


Shimada Seamount 11) 
Society Ridge [9] 


St. Helena [31] 
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Components 


Loihi [15] 
Mauri [3] 
Molokai [5] 
Oahu [17] 


Heard Island [9] 
Kerguelen Island [20] 
Kerguelen Plateau [12] 


Mauritius [1] 
Réunion [7] 


Manu a [4] 
Savai'i [8] 
Tutuila [9] 
Upolu 113) 


san Felix [4] 
san Ambrosio [1] 


Mehetia [2] 
Moua Pihaa [1] | 
Tahaa [1] 
Teahıtıa 14) 
dredge 11) 


Keferences 


27 
125005 
20 O 

6 

O 

IES 
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Table 1.1. Continued. 


Feature Components References 
Trinidade [1] ] 
Tristan de Cunha [5] 1522 
Tubuai-Austral Islands [22] 
Marotiri [1] 2 
Raevavae [1] ] 
Rapa [3] 52: 
Rimatara [4] | 
Rurutu [4] 21,25 
Tubuai [9] 5 
Walvis Ridge [10] 24 


!In the reference column, a "?" indicates a sample with an unknown reference. 


Reference guide: [1] Allčgre et al., 1987; [2] Barling and Goldstein, 1990; [3] 
Castillo et al., 1988; [4] Chaffey et al., 1989; [5] Chauvel et al., 1991; [6] Cheng 
et al., 1988; [7] Cohen and O'Nions, 1982a; [8] Davies et al., 1989; [9] Devey et 
al., 1990; [10] Duncan et al., 1986; [11] Dupuy et al., 1987; [12] Gautier el al., 
1990; [13] Gerlach et al., 1987; [14] Gerlach et al., 1988: [15] Gerlach et a/., 
1986; [16] Graham, 1987; [17] Halliday et al., 1990; [18] Halliday er al., 1988; 
[19] Hart, 1988; [20] Hart, unpublished; [21] Nakamura and Tatsumoto, 1988; 
[22] Newsom et al., 1986; [23] Palacz and Saunders, 1986; [24] Richardson e! al., 
1982; [25] Roden, 1982; [26] Salters, 1989; [27] Staudigel et al., 1984; [28] Stille 
et al., 1986; [29] Stille et al., 1983; [30] Storey et al., 1988; [31] Taras and Hart, 
1987; [32] Tatsumoto, 1978; [33] Vidal et al., 1984; [34] Weis, 1983; [35] Weis et 
al., 1987; [36] Weis et al., 1989; [37] West et al., 1987; [38] White, unpublished; 
[39] White and Hofmann, 1982; [40] White et al., 1989; [41] Dosis: ad and 
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Table 1.2. Isotopic means and standard deviations! for the OIB and the 


OIB+MORB data sets. 


Sr 
0182 
Меап 0.703943 
Std Dev 0.000892. 
OIB+MORB? 
Mean 0.703752 
Std Dev 0.000936 
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Nd 


0.512825 
0.000145 


0.512869 
0.000170 


6/4 Pb 


19.065 
0.880 


18.939 
0.870 


İlsotopic variance is the square of the standard deviation. 


?Mean and standard deviation based on 477 samples. 
Mean and standard deviation based on 567 samples. 


7/4Pb 


Де 
0.093 


1397] 
0.093 


8/4Pb 


38.965 
0:093 


355/99 
0.748 
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CHAPTER 2 


MATHEMATICAL AND STATISTICAL METHODS OF DATA ANALYSIS 


INTRODUCTION 

When dealing with a multidimentional data set with dimension greater 
than three, it is impossible to visualize the shape of the data in that space. This 
makes it difficult to choose "end-members" for the data, where end-members are 
interpreted as the vertices of the smallest simplex, with linear or nonlinear edges, 
that completely encloses all the data points. Previous work using two- 
dimensional plots to estimate the groups or end-members (Zindler et al., 1982; 
White, 1985; Zindler and Hart, 1986) can be misleading since those plots are 
projections of a higher-dimensional shape. For this study, it is possible to reduce 
the dimensionality of the OIB+MORB data set, via principal component analysis, 
and still retain its general shape, making it possible to choose end-members in 
three-dimensions. 

For the OIB data set, the data locations [oceanic islands] are not distributed 
evenly about the globe. This prompts the question as to whether there is any 
relationship between location and isotopic signature. To address this, a spatial 
correlation test (Mantel, 1967) is used to test for a correlation between the 
geographic distance and the "isotopic distance" between samples. In addition, a 
count is kept of the number of times a sample's isotopic "nearest-neighbor" 
occurs within the same island and within the same geographic feature. 

Finally, the globe has been divided by Hart (1984, 1988) into the islands 
lying inside the DUPAL belt, from 2? S to 60? S, and those lying outside. To see 
if there is statistical justification for separating the data into these two different 
populations, isotopic nearest-neighbor discriminant analysis 1s performed on the 


data set to obtain a misclassification error rate. The significance of this error 
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rate is based upon a randomization test of Solow (1990). While giving 

promising results, the randomization test for significance 1s inconclusive because 
spatial correlation within geographic features has not been accounted for. As an 
alternative, discrimination between isotopic signatures on the scale of geographic 


features inside and outside the DUPAL belt is addressed graphically. 


PRINCIPAL COMPONENT ANALYSIS 
Theory 
Principal component analysis can be viewed as a coordinate system 
transformation, but one that has particular properties. It generates a new set of 
variables, the principal components, that are linear combinations of the original 
variables: 


Zi = ES Es Xj — RES 


where the Z;s are the principal components, the ejjs are the transformation 
coefficients, and the Wë are the original isotope measurements (X, — 87 sr/86gr, 
Xor 143Nd/144Nd, X= 206Pp/204Pp, IE 207pp/204pp. Xs = 208pp/204Pb). 
The principal components have the following properties: 


(1) Z; and Z; are uncorrelated, for all 1, / 


(2) Variance(Z,) 2 Variance(Z5) 2... 2 Variance(Z5) 


3 
. 2 
(3) for all i, >» čje ] 
yel 


The transformation coefficients are the elements of the unit eigenvectors of the 5 
x 5 data covariance matrix. Because the isotopic ratios are on different scales, 


the data set must be normalized in order for all of the isotopes to be treated 
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equally in the analysis. One way to do this is to take each sample and for every 


isotope subtract the mean and divide by the standard deviation (Table 1.2): 


where Xij is the jth isotopic ratio for the ith sample, etc. This method weights 
the information provided by all five isotopes equally. Alternatively, Allègre et 
al. (1987) develop their own empirical norm, the "gcologic norm", that takes 
analytical errors into account and is designed to give equal weight to all isotopes 


except 207pp/204pb, which has the largest analytical error. 


Application to the OIB+MORB Data Set 
Because DMM [depleted MORB mantle] is onc of the proposed mantle 


end-member components, I have chosen to do principal component analysis using 


all of the oceanic island data [477 samples] plus a wide selection of MORB data 


[90 samples]. The covariance matrix for the OIB--MORB data set and its 
eigenvectors and eigenvalues are shown in Table 2.1. The sum of thc 
eigenvalues is the trace of the covariance matrix, ie. the sum of thc diagonal 
elements. This is equal to 5 because the diagonal elements of the covariance 
matrix, the scaled isotope variances, are all 1. To find out how much of the 
variance of the scaled data set is accounted for by each eigenvector, and thus each 
principal component, divide the corresponding eigenvalue by 5. The first three 
principal components account for 97.596 of the variance of the data set. 
Therefore it 1s reasonable to use the three-dimension principal component data 
set to select end-member components. This has important implications for the 


OIB+MORB data set. In n-dimensional space, the polygon containing the fewest 
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vertices [7+1] is a simplex. Thus, the OIB+MORB data set would require six 
end-member components to completely define it, if it spanned the entire five- 
dimensional space. The fact that it can be adequately represented in three- 
dimensions implies that the OIB+MORB data set requires only four end- member 
components. 

A comparison of eigenvalues and corresponding percentages of variance 
from this study and from Allégre et al. (1987) for OIB+ MORB and OIB data 
sets is presented in Table 2.2. It should be noted that the OIB eigenvalues from 
this study are found using a Separate covariance matrix derived from the 477 
OIB samples alone, as is done by Allegre et al. (1987). Their analysis yielded 
similar results for a three-dimensional fit to the data [OIB+MORB: 99.2% versus 
97.595; OIB: 98.896 versus 97.3%]. Part of the small difference that does exist 
may be due to the fact that they used a smaller data set [OIB+MORB: 91 samples 
versus 567 samples; OIB: 53 samples versus 477], in addition to the different 
methods used to scale the data. 

The procedure outlined above for computing principal components is 


compacted into matrix form, Z = EY, with exact solutions: 


Z 
ОД Гепепевенев 1 wl 
ДИ. ә Zay İml : U allen. Yan 
EE ZAN кәс ә EE PP YAN 
251 Y УД, UA | na 7 


where N = the number of samples [567], the Y;;s are the normalized isotopic 


values and the eigenvectors are the rows of the matrix E. Three two- 
dimensional plots of the first three principal components, with general end- 
member regions indicated, (Figs. 2.1, 2.2 and 2.3) are presented for comparison 


with those of Allégre et al. (1987) (Fig. 2.4). Plotting the principal component 
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values for the samples versus each other is the same as plotting the projection of 
the OIB--MORB population onto its eigenvector planes as they have done. “The 
two sets of plots are very similar, but mirror images of each other. This is 
simply because the eigenvectors used were of opposite sign, in no way affecting 


the validity of either set of plots. 


Mantle End-Member Components 

In three-dimensional space the principal component data form a 
tetrahedron (Fig. 2.5). It should be noted that the tetrahedron is not aligned with 
the principal component axes, so two-dimensional plots of the principal 
component data do not give an exact indication of the location of the extreme 
points. End-member component values are chosen by eye at the extremes of the 
tetrahedron using a rotating three-dimensional plotting program. 

First, the "nonlinear" end-member points are chosen, those that just form 
the vertices of the tetrahedron (Table 2.3). These end-members are referred to 
as "nonlinear" because they define the vertices of the smallest simplex enclosing 
the data points which has linear and nonlinear edges. In geometry, a simplex ts 
defined as a polygon with planar faces, but I am extending this definition to 
encompass polygons containing nonplanar faces as well. The purpose of 
choosing particular end-member points 1s to be able to express all of the sample 
points as a combination of the four end-member components, for later usc in 
spherical harmonic expansions. Though linear mixing 1s believed to exist 
between HIMU and EMI (Hart et al., 1986) and HIMU and DMM (Hart. 1988), 
more complicated mixing arrays are probable amongst the other components. 
Since no models exist for the nonlinear mixing arrays, it is easiest to represent 
the sample points as a linear combination of the end-member points. Thus, it is 


necessary to find the vertices of the smallest simplex with planar faces that 





28 


encloses as many data points as possible; these vertices are the "linear" end- 
members. These end-members are chosen by rotating the figure to look at the 
four sides of the tetrahedron edge on and moving out the "nonlinear" end- 
members until the planar-sided tetrahedron defined by linear mixing expands to 
contain as many sample points as possible, without becoming overly extreme 
(Table 2.4). This 1s an admittedly subjective process, but more accurate than 
choosing end-members using two dimensional plots. Figures 2.6 - 2.9 show the 
four views normal to each of the tetrahedron faces. 

When assuming linear mixing, the simplex defined by the final chosen 
"linear" end-member points excludes only 13 OIB data points, out of 477, and 3 
MORB data points out of 90 (Table 2.5), compared to the 85 OIB and 49 MORB 
data points excluded when using the “nonlinear” end-member values. The 
excluded points will have negative amounts of some of the end-members and will 
not be used in spherical harmonic expansions. 

The end-member values selected in principal component space are 
converted back into normalized isotope values (Tables 2.3 and 2.4) by 
substituting zeros [the mean value for each principal component] for the fourth 


and fifth principal components in the Z matrix: 


C DEER C14 ES Zu 
e11 €12€13€14€15 Hİ 

E eg Cr 20350. Z24 

Сб лил. Ca izi : i U Z34 | 

C41 5 C44 €Es1 €s2 €s3 €54 €55 155. O | 

бак: Сәд 025. 0 J 


where C; is the normalized 87Sr/80Sr ratio for the ith end-member component, 


and so forth. There is some error involved in this process, but because the 
variances of the fourth and fifth principal components are small, the error is 


small. To compute these errors, the entire OIB+MORB data set is transformed 
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into principal components; the fourth and fifth principal components are 
dropped; the approximate normalized isotope values are computed as above; and 
these values are then unnormalized and compared to the actual isotope values. 
The average absolute errors for this transformation are fairly small compared to 
the isotope standard deviations (Table 2.6). Compared to the range of analytical 
errors, all of the transformation errors are reasonable except the one for 
206pp/204pp, vvhich is approximately 6x larger than its analytical error (Table 
7) 

Finally, the samples are computed as percentages of the four “linear” end- 


members: 


1105 C14 Yy 
(51 Алы, Cr Pij Y); 
Ae Cry Paj| — Y3j 
Caras C44 P3j Y 4j 
(ас иә Си Рај Y 5j 
T EE 1d La 


where Pij is the percentage of the ith end-member component for the jth sample 


and Yi; is the ith normalized isotope value for the jth sample. The C matrix is 


the normalized end-member isotope value matrix computed from above with an 
additional row of ones. This row of ones and the one included in the Y vector 
define a constraint that the sum of the percentages add up to 1. This is necessary 
to provide useful positive results between O and 1 since the tetrahedron is not a 
four-component composition diagram, but resides in Euclidean space. QR 
decomposition 1s used to solve this over-determined system of equations. İt 
decomposes the C matrix into two matrices: Q [orthogonal] and R [upper 


triangular]: QRp= Y, with solutions: p*S! = R-İQİY. 
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SPATIAL CORRELATION TESTING 
Methodology 
In order to check for spatial correlation, a paired distance approach is 
employed, as outlined in Mantel (1967), using geographic and isotopic distances. 
The geographic distance used is that of an arc on a sphere connecting any two 


sample locations, ie. a great circle distance (Turcotte and Schubert, 1982). The 


angle Aij between the two locations / and J on the sphere (Fig. 2.10) is given by: 
Aj; = cos:! [cos 0j cos 0; + sin O; sin O; cos Kenn 


where 0; and q; are the colatitude and longitude of location / and 0; and q; are 
the colatitude and longitude of location J. The surface distance s between / and J 
IS: 

Sij — RA; 


where œR is the radius of the earth [R = 6378.139 km]. The isotopic distance used 
is the generalized Euclidean distance in multidimensions scaled by the variances 
of the isotopic ratios. Scaling by the variances of the isotopic ratios is necessary 


to keep the distance measurement from being dominated by the isotopic ratio 


with the largest variance, 206pp/204 pp (Table 1.2). For any two samples X; and 


X j the isotopic distance between them, d, is: 
dr s= Y (X;- X)I V (X;- X; 
1) (A; j) ( 1 j) 


where 
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is the isotope vector for ith sample [X1; 1s the 87Sr/86Sr ratio of the ith sample, 


etc.] and V is the diagonal variance matrix. A similar distance measurement, 
called Mahalanobis distance (Manly, 1986) was considered, but not used because 
it utilizes the covariance matrix. Covariance is a meaningful measurement when 
the data is normally distributed (elliptical) in space. From the three-dimensional 
principal component plots (Figs. 2.5-2.9), it is apparent that the data set is not 
elliptical, so covariance is a meaningless measurement concerning the nature of 
the data. 

Next, the correlation between the two distances for all the samples is 
calculated. The key to Mantel's (1967) technique is to determine the significance 
‘of this observed correlation by creating random pairings of the sample locations 
and isotopic signatures, calculating the appropriate distances, and computing 
their correlation, thus constructing a distribution against which the observed 
value can be judged. This distribution is that of the correlation under the null 
hypothesis that the geographic distances are matched to the isotopic distances at 
random. 

Zindler and Hart (1986) noted a relationship between the scale length of a 
geographic feature and the isotopic range of that feature. Basically, they 
concluded that the largest isotopic ranges occur in the largest geographic 
features, while small isotopic ranges may occur in small or large features. This 
implies a correlation between the within-feature geographic distance and the 
within-feature isotopic distance. The paired distance correlation method outlined 


above computes the correlation between geographic and isotopic distances both 
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within features and between features. In using this method, it 1s possible that any 
correlation within the features may be masked by a lack of correlation between 
the features. As an additonal check for within-feature correlation, a count 1S 
kept of the number of times a sample's isotopic nearest-neighbor [the sample that 
is the smallest isotopic distance from the sample in question] occurs within the 
same island and within the same island group [or island, if an island is not part of 
a larger group]. The counts are performed both for the observed data and for 
the random permutations. Those from the random permutations can be used, as 
before, to judge the significance of the observed counts. The larger scale 
geographic divisions of the data set into island groups and the remaining solitary 


islands (Table 2.7) will be referred to from this point on as features. 


Application to the OIB Data Set 

For this application, the OIB data is used since only oceanic island 
interrelationships are of interest. Two 477 x 477 distance matrices are calculated 
for the geographic and isotopic distances between samples. For the observed 
data, the correlation between the distance matrices is 0.1756 and the within island 
and feature nearest-neighbor occurrence rates are 61.4% and 76.7%, 
respectively (Table 2.8). The occurrence rates within islands and features appear 
significant and are confirmed so by randomization, as none of the generated 
occurrence rates are as large as the observed rates for 100 permutations (Table 
2.5). The correlation, on the other hand, is small, but attains significance 
compared to the randomization values which are all less than the observed value 
(Table 2.8). Thus, both methods indicate that there is spatial correlation between 
sample location and isotopic signature and the correlation that exists between 


samples within the same geographic feature seems to dominate. 
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Treating the samples inside and outside the DUPAL belt separately and 
then testing for spatial correlation yields results similar to those obtained with 
the whole data set (Table 2.8). 

It is not clear if all of the spatial correlation is due to the correlation 
within the features. There may be some additional spatial correlation between 
features. To check this, the appropriate samples are averaged to get an average 
isotopic signature and location for each feature (Table 2.7). Using all of the 
features both inside and outside the DUPAL belt, the observed correlation 1s 
0.1584 with a significance level of 0.13 [there are 13 permutations, out of 100, 
that have correlations higher than the observed correlation] (Table 2.8). Thus, it 
appears that there is spatial correlation between features. However, if there is a 
distinction between features inside and outside the DUPAL belt, this distinction 
may manifest itself as spatial correlation when testing all of the features at once. 
Testing the features inside and outside the DUPAL belt separately results in 
correlations of 0.0685 and 0.2645 with significance levels of 0.95 and 0.51, 
respectively (Table 2.8). These results indicate that there is no significant spatial 
correlation between the features, but that there is a distinction between features 


inside and outside the DUPAL belt. 


DISCRIMINANT ANALYSIS 
Isotopic Nearest-Neighbor Discriminant Analysis 
Methodology. Without taking account of spatial correlation, the validity 
of the division of the OIB data into samples inside and outside the DUPAL belt is 
addressed using isotopic nearest-neighbor as a discrimination rule. Using the 
isotopic distance measure outlined earlier, a given sample's isotopic nearest- 


neighbor is the sample that is the smallest isotopic distance away. 
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For the discriminant analysis, the assumption is made that the selected 
sample's location is unknown, so it is assigned the location of its isotopic ncarest- 
neighbor. This assigned location is compared to the actual location; if they are 
different, it is a misclassification. A count is kept of the number of 
misclassifications to calculate an error rate. 

Solow (1990) proposes a randomization technigue for judging the 
estimated misclassification probability or error rate. The importance of the 
misclassification error rate is to test the null hypothesis that there 1s no 
difference between the samples inside and outside the DUPAL belt. This is a 
trivial matter if the sampling distribution of the error rate under the null 
hypothesis is known, but in this case it is not. A simple but effective way to 
judge the significance of the observed error rate is to construct a randomization 
distribution under the null hypothesis that the pairing of isotopic signatures and 
locations inside or outside the DUPAL belt occurs by chance. Applying the 
randomization technique to the data, the samples retain their isotopic signature, 
so their isotopic nearest-neighbor remains the same, but they are randomly 
assigned to locations inside and outside the belt. The discriminant analysis is 
done, as described above, with this new randomly constructed data set to get its 
misclassification error rate. Then the process is repeated to construct the 
distribution. 

Application to the OIB data set. For the OIB data set, the observed 
misclassification rate is 7.3% and the randomization error rate ranges from 
35.2% to 53.7%. Superficially, it appears that describing the data as two 
populations residing inside and outside the DUPAL. belt is viable. However, the 
within-feature spatial correlation has not been accounted for in this analysis, If 
16.1% of the time, a sample's isotopic nearest-neighbor is located within the 


same geographic feature, then it seems obvious that the misclassification error 
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rate would be small. The observed error rate itself is not incorrect, but the 
randomization distribution of error rates against which it is being judged is 
incorrect. In order for the significance of the observed error rate to be properly 
judged, the spatial correlation must be preserved in the randomization process. 
In this case, preserving the spatial correlation is too complicated to pursue when 


other methods may provide the desired information. 


Graphical Discrimination of Geographic Features 

As shown earlier, the correlation between isotopic distance and the 
geographic distance within features is very strong. A way around this spatial 
correlation is to look for differences between populations inside and outside the 
DUPAL belt on the feature level. The averaged isotopic values for the features 
(Table 2.7) are scaled by the mean and standard deviation of the isotopes derived 
from the entire OIB+MORB data set (Table 1.2) and expressed in terms of 
principal components using the eigenvectors of the OIB+MORB correlation 
matrix (Table 2.1). 

The first three principal components are plotted to Jook for differences in 


features inside and outside the DUPAL belt, with the general direction of the 


end-members indicated (Figs. 2.11-2.13). On all of the plots, but especially Z4 


versus Z^, most of the features outside the belt cluster in a band between DMM 


and HIMU, with the exception of the Hawaiian islands [the Koolau volcanics on 
Oahu show a strong EMI signature (Hart, 1988)], Shimada Seamount [which also 
has an EMI signature (Hart, 1988)], and the Azores [Sáo Miguel has a strong 
EMII signature (Hart, 1988)]. Essentially, the features outside the DUPAL belt, 
with few exceptions, occupy only part of the available isotopic space, while 
features inside the belt occupy all of the available isotopic space, including some 


overlap with features outside. This is essentially the relationship found by Hart 
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(1988), not that the two populations are totally separated, but that one population 
contains isotopic signatures that the other does not. This unportant that this two 
population distinction is still valid on tlie feature level, Since itis sti apparent ` 
at this larger scale [not just sample to sample] the geochemical signatures of the 
oceanic island basalts do have a lonp, wavelength component to then, making it 
feasible to attempt to quantify these signatures using spherical harnionie 
expansions. 

In addition to this grapliical presentatiou, the diserinination analysis can 
also be done on the feature level, but the variances of the isotopes within each 


feature must be accounted for in some way. 


SUMMARY 
Mathematical and statistical methods to explore and characterize the OIB 


and MORD data reveal these main points: 


« OIB+MORB data require only four mantle cnd member components to 


completely span the range of known isotopte values. 


© Choosing the mantle end-tnember components can be made easier [and 


more accurate] with the usc of principal component analysis, 


. Within geographic features, there is a significant correlation between 
location and isotopic signature, but between geographic features, there 


Is not. 
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Graphical discrimination of geographic features shows that the 
distinction between islands inside and outside the DUPAL belt is 


viable. 


The existence of the DUPAL anomaly on the feature level indicates 


that the anomaly has a long wavelength component to it. 
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Table 2.1. Covariance matrix! of the five isotopes with its eigenvectors and 


eigenvalues. 


Isotope 


Y1 
Y2 
Y3 
Y4 
Ys 


Isotope 


Yi 


1.000000 


-0.017647 
-0.122196 
0.565079 
0.574974 
0.578661 


Eigenvalues 2.830 


Covariance Matrix 


Y2 


-0.796442 
1.000000 


Y3 


-0.273004 
0.054987 
1.000000 


Уд 


-0.019107 
-0.170370 
0.901205 
1.000000 


Eigenvector Matrix 


II 


0.699432 
-0.683457 
-0.195210 
-0.006763 

0.074352 


1.861 


MI 


0.682362 
0.679315 
0.019019 
0.249290 
-0.102014 


0.183 


IV 


-0.174033 
-0.179201 
-0.235936 

0.753098 
-0.561051 


0.09] 


Ys 


0.061599 
-0.295078 
0.894577 
0.901429 
1.000000 


-0.120738 
0.156120 
-0.765867 
0.200145 
0507 


0.035 


Percentage of total variance accounted for by each eigenvector 


56.6 


5 


EEN 


1.8 


0.7 


l Only the upper half of the covariance matrix 1s shown since it is symmetric. 


All eigenvector values are rounded to six decimal places from the fourteen 


decimal accuracy used in the calculations. 


Covariance matrix is calculated using 477 OIB and 90 MORB samples. 
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Table 2.2. Comparison of eigenvalues and Percentages of variance accounted for 
by the corresponding cigenvectors from this Study and from Allegre eral. 
(1987)! for OIB+MORB and OIB data sets. 





OIB+MORB I II III IV V 

2 2.830 1.861 0.183 0.09] 0.035 
[56.6%] [37.2%] [3.7% | [1.8%] [0.794 

2 Allégre 3.20 1.6] 0.15 0.03 0.01 

et al. [64.095] [32.2% | [3.0% | [0.6%] 10.2%] 

OIB I II HI IV V 

4 3.047 1.568 0.249 0.099 0.037 
[60.9% | [31.4%] [5.097 12.0%) 10.7%] 

Allegre 2.85 1.87 (12: 0.05 0.01 

et al. [57.0%] [37.4%] [4.4%] [1.0%] 10.2%] 


 —rñÑsas ss nn 


lEigenvalues from Allégre el al. (1987) are converted to scaled eigenvalues that 
add up to 5 for comparison with eigenvalues from this study. 


Percentages of variance accounted for by the corresponding eigenvectors are 
indicated in parentheses. 


ŽBased on 367 samples. 
Based on 91 samples. 
“Based on 477 samples. 
Based on 33 samples. 


I —_AAAA= “ r.o - ——— —GÀ— — 
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Table 2.3. "Nonlinear" end-member component values in principal component 


space and the transformed values in isotope space. 


End-Members in Principal Component Space 


Z] 72 zn 
EMI 7:7: 3:6 1.3 
EMII 1.0 4.0 27 
HIMU 5.0 Ee >> 
DMM 77 2% 0.4 


— um e m 


End-Members in Isotope Space 


EMI 
EMI 
HIMU 
DMM 


XI 


0.705311 
O TOVO 
0.702659 
0.702171 


X2 


0.512343 
0.512638 
0.512887 
0.513329 


X3 


1 
18.788 
21615 
17.594 


X4 


15.43] 
15.07 
1200 
15.38] 


X 5 


38.232 
207287 
40.911 
36.983 
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Table 2.4. "Linear" end-member component values, based upon linear mixing, 
in principal component space and the transformed values in isotope space. 


End-Members in Principal Component Space 


Z1 72 73 
EMI -2.4 3.6 -1.6 
EMI 1.8 4.5 2.6 
HIMU 6.0 -1.9 -0.6 
DMM -4.3 -3.7 5 


End-Members in Isotope Space 


X1 X2 X3 X4 X 5 
EMI 0.705126 0.512316 17.121 15.403 38.082 
EMII 0.708329 . 0.512609 19.103 15.724 27650 
HIMU 0.702026 0.512896 223203 155 41.337 
DMM 0.701624 0.513428 17.459 15.351 36.704 
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Table 2.5. Samples excluded from linear mixing tetrahedral volume.! 





Location Sample Number Row Number? 
Azores, São Miguel SMID Do 
SM49 36 
Galapagos E55 173 
Gough 10 175 
Hawaii 69Tan2 200 
Kerguelen Plateau DR02/12 27 
DROS 282 
747c-12r-4-45-46 27 
747c-16r-2-81-84 294 
Marquesas uapll 329 
Pitcairn, Pul wana 642 370 
St. Helena 7 469 
237 482 
Atlantic Ocean AD3-3 535 
SVV Indian Ridge DI 536 
D5 559 


lOIB samples excluded from the volume will not be used in spherical harmonic 


expansions. 


?Indicates row number of the data set included in Appendix A. 
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Table 2.6. Average absolute errors in transforming three-dimensional principal 
component data into five-dimensional isotope data, with their ratio to isotope 
standard deviations and comparison to analytical errors. 


X1 X2 X3 X4 A 


/Ј1 


Average absolute error [for 567 samples] 
0.000041 0.000008 0.111 0.016 0.111 


Ratio of average absolute error to standard deviation! 


0.043802 . 0.047124 0.128 0.170 0.149 
Absolute Error Percentage Range? Analytical Error Range 
X1 0.00580 to 0.00584% 0.003 to 0.01% 
X2 0.001558 to 0.001561% 0.003 to 0.01% 
X3 0.225 to 0.330%/amu 0.03 to 0.05%/amu 
X4 0.0337 to 0.0347%/amu 0.03 to 0.05%/amu 
X5 0.068 to 0.075%/amu 0.03 to 0.05%/amu 


İlsotopic standard deviations for the OIB+MORB data set are indicated in Table 
ри. 


2Absolute error percentage ran ges are calculated using the average absolute 
errors and the ranges of the isotopes in the OIB+MORB data set: 


Xl 0.702290 to 0.707400 
E92 0.512376 to 0.513290 
X3 16.943 to 21.755 
X4 15.406 to 15.862 
X5 37.235 to 40.619 
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Table 2.8. Correlations between geographic and isotopic distance matrices and 
island/feature isotopic nearest-neighbor occurrence rates for all samples in the 
OIB data set and samples inside and outside the DUPAL belt. Correlations for all 
the geographic features and those inside and outside the DUPAL belt are also 
given. 


Correlation Randomization Rangel 
OTD 0.1756 -0.0077 to 0.0224 
Inside DUPAL 0.0641 -0.0035 to 0.0343 
Outside DUPAL 0.6142 -0.0045 to 0.0699 


Island Occurrence Rate Randomization Range! 


OIB 61.4% 1.0% to 3.6% 
Inside DUPAL 63.6% 1.6% to 6.2% 
Outside DUPAL 47.1% 1.7% to 11.0% 


Feature Occurrence Rate Randomization Range! 


OIB 16.1% 3.496 to 10.096 
Inside DUPAL 71.590 4.2% to 11.1% 
Outside DUPAL 15.0% 1020006 


Correlation Randomization Range! Significance 


Level 
Features 0.1584 0.0369 to 0.2075 015 
Inside DUPAL 0.0685 0.0350 to 0.2807 0.95 
Outside DUPAL 0.2645 0.0649 to 0.6602 0.51 


İRandomization ranges based upon 100 random permutations. 








50 


i 
3 
1 
1 
3 
j 
1 
3 
1 
1 
1 
1 


22 





Fig. 2.1. Plot of the second principal component versus the first principal 
component for the OIB+MORB data set. Symbols: x = MORB data, open circle 
= OIB samples inside the DUPAL belt, black diamond = samples outside the 
DUPAL belt. General mantle end-member component regions are indicated. 
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Fig. 2.2. Plot of the third principal component versus the first principal 
component for the OIB+MORB data set. Symbols: x  MORB data, open circle 
= OIB samples inside the DUPAL belt, black diamond = samples outside the 
DUPAL belt. General mantle end-member component regions are indicated. 
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Fig. 2.3. Plot of the third principal component versus the second principal 
component for the OIB+MORB data set. Symbols: x  MORB data, open circle 
= OIB samples inside the DUPAL belt, black diamond = samples outside the 
DUPAL belt. General mantle end-member component regions are indicated. 
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Fig. 2.4. Plots of the second principal [V2] component versus the first [V1] and 
the third principal component [V3] versus the first [V1] for a smaller 
OIB+MORB data set from an analysis done by Allégre et al. (1987). These plots 
are the mirror images of the ones done for this analysis because the chosen 
eigenvectors for the two analyses are of opposite sign. 
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Fig. 2.5. Three-dimensional view of the OIB+MORB principal component data. 
Axes: Х = 21, У = Z2, Z = Z3. Symbols for the end-member components: += 
EMI, x = EMII, diamond = HIMU, square = DMM. 
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Fig. 2.6. Three-dimensional view of the OIB+MORB principal component data 
parallel to the EMI-EMII-HIMU plane. Symbols for the end-member 
components: + = EMI, x = EMII, diamond = HIMU, square = DMM. 


ə 
» Y 
b ' 
E ° 
b 
, 
P 
( 
4 
: 1 D pari 
? ma?! fie? Ose 
A" 1 1 p 





sle 





Fig. 2.7. Three-dimensional view of the OIB+MORB principal component data 
parallel to the EMI-EMII-DMM plane. Symbols for the end-member 
components: + = EMI, x = EMI], diamond = HIMU, square = DMM. 
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Fig. 2.8. Three-dimensional view of the OIB+MORB principal component data 
parallel to the EMI-HIMU-DMM plane. Symbols for the end-member 
components: + = EMI, x = EMII, diamond = HIMU, square = DMM. 








Fig. 2.9. Three-dimensional view of the OIB--MORB principal component data 
parallel to the EMII-HIMU-DMM planc. Symbols for the end-member 
components: + = EMI, x = EMII, diamond = HIMU, square = DMM. 
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Fig. 2.10. Geornetry for determining the surface distance s between locations I 
and J on the globe, where 6 and q are colatitude and longitude and A is the angle 
between the two locations taken from the center of the Earth, From Turcotte 
and Schubert (1922). 
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Fig. 2.11. Plot of the second principal component versus the first principal 
component for the 36 geographic features. Symbols: open circle < features 
inside the DUPAL belt, black diamond = features outside the DUPAL belt. 
Labeled points: 1 = Hawaiian Islands, 2 = Shimada Scamount, 3 = Azores. The 
general directions of the mantle end-member component regions are indicated. 
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Fig. 2.12. Plot of the third principal component versus the first principal 
component for the 36 geographic features. Symbols: open circle — features 
inside the DUPAL belt, black diamond features outside the DUPAL belt. 
Labeled points: 1 = Hawaiian Islands, 2 = Shimada Seamount, 3 = Azores. The 
general directions of the mantle end-member component regions are indicated. 
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Fig. 2.13. Plot of the third principal component versus the second principal 
component for the 36 geographic features. Symbols: open circle = features 
inside the DUPAL belt, black diamond = features outside the DUPAL belt. 
Labeled points: 1 = Hawaiian Islands, 2 = Shimada Seamount, 3 = Azores. The 
general directions of the mantle end-member component regions are indicated. 
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CHAPTER 3 


SPHERICAL HARMONIC REPRESENTATION OF ISOTOPIC 
SIGNATURES: THE CONTINUOUS LAYER MODEL 


INTRODUCTION 

Hart (1984) contoured world maps of OIB isotope data for his three 
DUPAL anomaly criteria (ASr » 40; A7/A » 3; A8/A » 40]. These maps show a 
concentrated band spanning approximately 60? of latitude, centered on 309-4095, 
with pronounced highs for the anomaly criteria in a region from the South 
Atlantic to the Indian Ocean [ASr, A7/4, A8/4] and in the central Pacific [ASr, 
A8/4]. Qualitatively, Hart (1984, 1988) believes this geochemical anomaly 
correlates with other geophysical anomalies: the slab-corrected geoid (Hager, 
1984), deep mantle P-wave tomography maps (Dziewonski, 1984), slow P-wave 
regions at the core/mantle boundary (Creager and Jordan, 1986) and equatorial 
anomalies in the core (Le Mouél et al., 1985). These geophysical anomaly 
patterns are typically expanded in terms of spherical harmonics, therefore any 
attempt to make a quantitative comparison between geochemical and geophysical 
patterns requires expanding the geochemistry data in spherical harmonics as 
well. 

Expansion of the geochemistry data is approached in two ways, based 
upon an assumed geometry for the OIB geochemical reservoir.: The first 
approach, the "continuous layer model" discussed in this chapter, assumes that 
the OIB reservoir is a continuous layer [not ruling out heterogeneities within this 
layer] and tries to reconstruct this layer. Plumes from this layer only sample the 
continuous geochemical "function" in discrete locations. With the geochemistry 
"function" unknown, the spherical harmonic coefficients must be solved for 
using least squares, singular value decomposition or a similar method that will 


approximate the values of the geochemistry "function" where there is no data. 
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The second approach, the "delta-function model" discussed in Chapter 4, assumes 
that the OIB reservoir is composed of a series of point sources, each feeding a 
separate plume. In this case, the geochemistry "function" is known and can be = 
represented as a series of delta-functions. The spherical harmonic coefficients 
can be solved for directly with the simplification from integration to summation 
allowed by the delta-function approximation. 

The continuous layer model and the delta-function model are not meant to 
suggest two end-member possibilities for OIB source geometry. Rather, the 
delta-function model can be regarded as an approximation of the continuous 
layer model that gives a mathematically robust solution for the spherical 
harmonic coefficients. In regard to the oceanic crust model of Hofmann and 
White (1982), the continuous layer model corresponds to the accumulated layer 
of subducted oceanic crust, with the plume-forming instabilities occurring at 
discrete locations within this layer. The delta-function model can also be 
reconciled with the accumulated layer model, with the stipulation that discrete 
pockets [point sources] within this layer form and feed individual instabilities. 

For the purposes of minimizing small scale variations [1e. variations 
within a single island or island chain] in the geochemistry "function" that cannot 
be accurately represented with the incomplete global data coverage, this spherical 
harmonic study is based on the averaged isotopic signatures of the 36 geographic 
features (Table 2.6). These average isotopic values are converted to mantle-end 
member component percentages (Table 3.1), as outlined in Chapter 2, to form 


the data matrices used in the expansions. 


SPHERICAL HARMONIC BASICS 


Spherical harmonics, Y/'(6,@), are a set of orthonormal functions over the 


unit Sphere: 
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(2/+1)U-m)! 


a) P7 (cos) eim 


Y¡ (0,9) = 


where / is the degree of the expansion, » is the order of the expansion, 0 is 
colatitude [0 = xx/2 - lautude; O<0<n] and @ is longitude [-T<Y<T). The 

functions e! form a complete set of orthogonal functions in the index m on the 
interval -T<p<T and the associated Legendre polynomials P7'(cos0) form a 

similar set in the index / for each 7: value on the interval -1<cos0<1 (Jackson, 
1975). Therefore their product forms a complete orthogonal set on the surface 
of the unit sphere in the two indices /,m. The spherical harmonic functions used 
in this analysis are normalized by the square root term so that their integrated 
Square over the sphere is unity [in most geophysics applications, the functions are 


normalized so that the integrated square over the sphere is 4x]: 


2x i 


аф | a(cos0)YT (0,9)' Y7(0,9) 2 87/8, 


where the asterisk denotes complex conjugation. 


Any function f(0,9) can be expanded in spherical harmonics: 


L l 
f(0.)) A > C Y (OO) 


l=) m=-l 


where L is the maximum degree of the expansion and C;' are complex spherical 


harmonic coefficients. Written in a more explicit form, the equation becomes: 


[T DRM 
IDEE 3 P/'(cos0) [ A7'cosnip *- B'sinmq| 
mE 4ru lam)! 
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where A7 and B/F are real Spherical harmonic coefficients. When expanding a 


function from degrees O to L, the number of coefficients that need to be 


L 


2l - 1 
calculated is: Ce There are actually an additional [L+1]} coefficients 


involved, but for m = О, sinm@ = 0, so p? < 0. [It is important to realize that 
only having 36 features limits the possible spherical harmonic expansion to 


degree 5, in order to avoid a purely underdetermined problem. 


MANTLE END-MEMBER COMPONENTS 

When attempting to use inverse methods to solve for the harmonic 
coefficients of an unknown function, careful attention must be paid to the 
variation of the data as a function of distance to avoid the problem of aliasing. 
For a simple two-dimensional case, aliasing occurs if the sampling interval is 
longer than half the shortest wavelength of the function sampled, causing the 
sampled points to show a periodicity that does not exist in the original data, The 
minimum distance between any two geographic features in the OIB feature data 
set is 33.396 km, but the distance between features is not constant. Plots of data 
variation versus distance between data locations make it possible to select a 
mininum sampling distance based on the shortest distance required to get the 
maximum data variation. This minimum sampling distance then controls the 
minimum degree to which the data must be expanded in order to adequately 
represent the data in spherical harmonics without aliasing. The relationship 


between wavelength and degree is: 
pese Te. 
y (1+1) 
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where À is the wavelength [A = 2"(sampling distance)], R is the radius of the 
earth [R < 6378.139 km] and / is degree. Solving for degree in terms of 
wavelength: 
-I+ AJ ] кај2ид| 

: A 


2 


i= 


Variation-Distance Relationships 
For variation-distance relationships, the distance measure is the angle Aij 
[in degrees] from the center of the earth between any two locations / and / [see 


Chapter 2] and the variation measure is the absolute value of the difference 


between the mantle component percentages at those locations. The angle Aij can 


be transformed into a great circle distance in km by converting Ajj to radians and 


multiplying by the radius of the earth R. 

Plots of absolute difference versus angle for the four mantle components 
(Figs. 3.1-3.4) show the maximum variation in the components occurring on 
very short distance scales for the EMI and HIMU components and moderate 
distance scales for the EMII and DMM components. Based upon these plots, the 
minimum sampling distances [in degrees] are ~ 14.5? for EMI and HIMU, ~ 39° 
for EMII and ~ 57° for DMM. These correspond to expansions out to degrees 
12, 4 and 3, respectively. For the current problem, the EMII and DMM data sets 
can be expanded in spherical harmonics as they are, but the EMI and HIMU data 


sets require some additional manipulation. 


Variation Reduction by Categorizing Features 
Separation of the geographic features into populations located inside and 
outside the DUPAL belt [2°S to 60°S] does not result in two distinct isotopic 


populations [Chapter 2]. Essentially, one population [outside the belt] defines a 





68 


small field in isotopic space, while the other population [inside the belt] defines a 
larger field that overlaps with the smaller field (Fig. 3.5). A possible source of 
the large, small-scale isotopic variation exhibited by the EMI and HIMU data sets: 
is the juxtaposition, due to the overlap in isotope space, of features having a 
strong DUPAL signature next to those that do not. If it is possible to separate 
DUPAL-type features [those features showing a strong DUPAL signature] from 
DMM-type features, this separation might reduce the small-scale variation within 
these two populations and thus reduce the degree to which the population data 
must be expanded. 

Since the goal is to separate DUPAL-type features from DMM-type 
features, a logical starting place is to look at the spatial distribution of different 
percentage categories of the DMM component in three-dimensional principal 
component space (Fig. 3.6). Six DMM percentage categories |<10%, 10-20%, 
20-30%, 30-40%, 40-50%, >50%] can be distinguished as six separate point 
groupings. Most striking is a large spatial separation that occurs within the 30- 
40% category for a small percentage difference [Louisville - 31.84%, Balleny - 
32.17%, Cocos - 38.53%]. This is a reasonable place to separate the DUPAL- 
type features from the DMM-type features, with a boundary value of 32% 
DMM, for simplicity. The resulting 27 DUPAL-type features and 9 DMM-type 
features, with their percentage of the DMM component are listed in Table 3.2. 

There are too few DMM-type features to draw any conclusions from plots 
of absolute difference versus angle. For the DUPAL-type features, plots of 
absolute difference versus angle of the DUPAL components [EMI, EMII and 
HIMU] show no reduction in the small-scale variation, while that of DMM does. 
with an increase in sampling distance from ~ 57° to 89° (Figs. 3.7-3.10). In 
retrospect, this is an obvious result of the artificial separation performed. The 


percentage categories are basically parallel slices through the tetrahedron that 
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move from a broad base of lower percentages to a peak of high percentages 
approaching an end-member component apex on the tetrahedron [like a ternary 
diagram]. It is true that these slices can separate DMM-type features from 
DUPAL-type features, but only the variation of the DMM percentages are 
reduced. To reduce the variation of the individual DUPAL components using 
this method, EMI-type features would have to be distinguished from non-EMI- 
type features, etc. This would generate four different, though overlapping, sets 
of features to use to characterize the four different components. Manipulation of 
the data set in this way is not desirable, so another method must be pursued in the 


attempt to reduce small-scale data variation. 


Variation Reduction by Filtering 

Another method to reduce small-scale variation [and hopefully enhance 
any long wavelength component] is to filter the data set in some way. Here, a 
simple circular filter, of fixed radius, is applied to each feature location. The 
new data values assigned to that feature location are the means of the mantle 
component percentages of the feature locations that fall within the circle. To 
ensure that there are always at least two features falling within the circle, the 
radius of this circle is determined by the longest distance to the nearest feature 
location. Nunivak Island is the most isolated feature with the nearest feature 
being the Hawaiian Islands at an angular distance [from the center of the earth] of 
40.86". The circle radius is then 40.9°, for simplicity. 

Plots of absolute difference versus angle for the filtered data set yield 
interesting results (Figs. 3.11-3.14). All of the mantle component data sets show 
a reduction in small-scale variation, except EMIT, which shows an increase in 
variation, with a decrease in angular sampling distance from ~ 39° to 27? 


[expansions to degrees 4 and 7, respectively]. The remaining plots show an 





70 


increase in angular sampling distance from — 14.5? to 37? [expansions to degrees 
Ipond respectively] for EMI, an increase from ~ 57° to 102° [expansions to 
degrees 3 and 2, respectively] for DMM and a dramatic increase from — 14.5? to 
83° [expansions to degrees 12 and 2, respectively] for HIMU. Now that the 
small-scale variation has been significantly reduced by filtering, the filtered EMI 


and HIMU data sets can also be expanded in spherical harmonics. 


INSIGHTS FROM GEOPHYSICAL DATA 

It is unclear how accurate the spherical harmonic expansions of the OID 
feature data set will be due to the limited global coverage and the highly variable 
nature of the data. In an attempt to address these problems, three geophysical 
data sets, with different variance characteristics, are constructed with the same 
limited coverage to provide a sort of control set against which qualitative 
comparisons can be made. Geoid, gravity and gravity gradient anomalies are the 
chosen geophysical measures because their coefficients are well known and they 
form a kind of continuum from the long wavelength [low degree] dominance in 
the geoid signature to the short wavelength [high degree] dominance in the 
gravity gradient signature (Fig. 3.15). Techniques applied to the mantle 
component data, to solve for the spherical harmonic coefficients, are also applied 
to these constructed data sets to see how closely the actual geophysical 


coefficients can be approximated. 


Construction of Geophysical Data Sets 
The gravitational potential V, in spherical harmonics as a function of 


radial distance r, is given by: 
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ии b 5.3 + 2. (R y » puo P" (cos0) | A cosmq + B fing] 


T 


where G is the gravitational constant [G — 6.6726х10- 11 m3/kges?], M is the 


mass of the earth [M = 5.973x1024 kg] and R is the radius of the earth in meters 
(Stacey, 1977). The gravitational potential anomaly [0V = Vobserved ` 


Viheoretical! IS: 


oul ett [CEDE | 
ÖV = Bux 2 k Sos! P7 (cos0)| ¡cosmo + Bj ə 


which can be converted to the geoid anomaly dN (in m) by dividing by 
2 =-GM/R2: 


ôN = 





= R Ver E e » pie P/ (cos0) | A/cosmqQ + “. 


T 


The geoid anomalies calculated here are referenced to a theoretical hydrostatic 
sphere to remove the effect of the earth's rotation (Hager, 1984). Gravity is the 
derivative of the gravitational potential with respect to radial distance, so the 


radial gravity anomaly is: 


2 ö(öV) 
BB or 
eZ | m m 
- š P pa R) 5 Fə P7 (cos0) | A"'cosmqQ + B! song] 
=? m= 


Gravity gradient is the derivative of gravity with respect to radial distance, so 


the radial gravity gradient anomaly is: 
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Ol = 9(Óg) 


dr 
co l 
BOM 1:7 m. EHEM)! PT (cos8) | AT cosmo + BY sinm | 
K > р2 i 2 dt lun a | ? d 





Evaluating at r = and using the spherical harmonic coefficients 2-20 from the 


GEM-L2 model (Lerch et al., 1982), the equations simplify to: 


20 1 
EXE ae E I 
EN = R (EINEM). pm (coso) lA See 
p 2 Ax (HR 7 (cos8) | ¡ cosmo + B; sinmq) J 


20 l 
Си ES > (1-1) > —. P; (cos) |A7 cosmo + - 
|=2 : 


m=() 


20 l 
77 OM > (I+1)(1+2) E (241) P; (cos) [A cosmo + m siamo)’ 
Ri İZ dml | 


m=() 
It is important to note that the GEM-L2 coefficients must be multiplied by ¥4x 
before they are plugged into these equations to be consistent with the spherical 
harmonic normalization used in this study. The three geophysical control data 
sets are constructed by calculating the values of the geoid, gravity and gravity 


gradient anomalies at the 36 feature locations (Table 3.3). 


Variation-Distance Relationships 

The different characteristics of the contructed geophysical data sets are 
apparent in plots of absolute difference versus angle (Figs. 3.16-3.18). The 
geoid plot shows a clean and fairly symmetric degree 2 pattern, with an angular 


sampling distance of ~ 102°. The gravity plot is a little more dispersed, with 
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weaker symmetry and an angular sampling distance of ~ 95° [expansion to degree 
2]. Finally, the gravity gradient plot shows even more dispersion and an angular 
sampling distance of ~ 67° [expansion to degree 3]. A comparison of these plots 
to those for the mantle components clearly illustrates the complexity of the 
geochemistry data. Even the gravity gradient data [dominated by short 
wavelength energy] appears to have less small-scale variation [larger angular 


sampling distance] than all of the mantle component data sets. 


Variation Reduction by Filtering 

The same circular filter technique outlined above is applied to the 
geophysics data to see its effect (Figs. 3.19-3.21). The filtered geoid data set 
retains its strong degree 2 signature [angular sampling distance ~ 93°], but there 
Is a slight increase in the dispersion of the data points. Like the geoid, the 
filtered gravity data maintains its angular sampling distance [~ 93°] and it shows 
a Slight decrease in data dispersion. The gravity gradient data is most affected by 
the filtering process. The data dispersion due to large variation at small and 
large angles is reduced. In addition, the angular sampling distance Is increased to 


~ 77°, corresponding to spherical harmonic expansion to degree 2. 


EXPANSION OF GEOPHYSICAL AND GEOCHEMICAL DATA SETS 

By choosing the sampling distances based upon the inherent variation- 
distance relationships of the different data sets, the problem of aliasing is 
eliminated. Of course, the location patterns that result from spherical harmonic 
expansions may not represent the true patterns as they exist in the mantle, but 
without a more extensive global data set, there is no way to better approximate 


the true pattern. Coefficients will be found for all six geophysical data sets 
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[filtered and unfiltered], for the EMIT and DMM data sets and for the filtered 
EMI and HIMU data sets. 

Solving for the spherical harmonic coefficients needed to expand a given 
function is a linear inverse problem. More specifically, the expansion of the 
mantle components or geophysical measures is a discrete linear inverse problem, 
since the data are discrete observations. The terminology and symbolog y used 
here to discuss inverse problems is that of Menke (1989). Values of the mantle 
components or geophysical measures at the feature locations form a vector of 
data values d [Nx1J. The unknown spherical harmonic coefficients form a 
vector of model parameters m [Mx1]. Relating the two is the data kernel matrix 
G [NxM], composed of Legendre polynomials [functions of colatitude] combined 
with sine and cosines [functions of longitude]. In matrix form the equation is: 


Gm = d, or written out more explicitly: 


y Po(cosBp) . P/(cos6p) cosLqQo VP; [(cos0o) sinLpo AO do 


2 | 
V Pi(cosö)y) . . . V Pi(cos6y) софи. Ҹ.Рӱ(созӦл/) sinLon x İN 
LL | 


where L is the maximum degree of the expansion, N is the number of data 


observations and Y is the normalization factor mentioned earlier. 


Least Squares Method 

Theory. If the equation Gm = d provides enough information to uniquely 
determine the model parameters or the best fit to the model parameters, then 
solving for the spherical harmonic coefficients from degrees O to 5 is an even- 


determined problem İN = 36, M = 36] and solving for the coefficients from 
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degrees O to <5 is an overdetermined problem [N < 36, M <36]. For an 
overdetermined system of equations Gm - d, with more equations than 
unknowns, there is no exact solution. The least squares method finds the model 


parameters that minimize the error between the observed data and the predicted 


data, ie. it minimizes the L5 norm of the prediction error: 


N 
[onoi | eb S » le,l” , Where e; = iR de 


i=] 

When solving for the model parameters m [spherical harmonic coefficients], it is 
best to use QR decompositon. The normal equations GİGm z G İd lead to the 
solution: m€$t - (GİG) İG Ta, but iif GİG is ill-conditioned, then taking its 
inverse leads to inaccurate solutions. QR decomposition 1s more accurate than 
the normal equations for ill-conditioned matrices. It decomposes the data kernel 
matrix G into two matrices: Q [orthogonal] and R [upper triangular]: QRm = 
d, with solutions: mêst = Hold 

Application. As a test of the viability of the least squares method, the 
spherical harmonic coefficients for the EMII percentage data and the geoid 
anomaly data are solved for in nested groupings from degrees 0-1 up to degrees 
0-5. As the data is expanded out to greater degrees, the coefficients should 
decrease smoothly. Table 3.4 shows how the degree 2 coefficients vary as the 
two data sets are expanded out to progressively higher degrees. Only the 2 and 
AS coefficients for the geoid and the Ağ coefficient for EMII decrease smoothly 
for the degrees 0-2 through degrees 0-4 expansions. The other coefficients 
either get larger or oscillate. When solving the even-determined system [degrees 
0-5], all of the coefficients experience a large increase or decrease, indicating a 


very unstable solution. 
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Since the geoid coefficients are known, the correlations [by degree] 
between the actual coefficients and the computed coefficients for the nested 
groupings can be calculated. The correlation coefficient rj for two sets of 
coefficients [A1,B1] and [A2,B2] is given by the ratio of covariance to variance 


at each harmonic degree (Richards and Hager, 1958): 


l 
S j o 02091713277] 
me m=U 
po l l 
SS SCTE (Kor 
m=0 m=0 


Correlations with the actual geoid coefficients can only be made at degrees 2 and 
higher since the actual degree O and | coefficients are zero. Correlations of the 


actual geoid coefficients to those calculated using least squares are: 


Expansion Correlation Coefficient. [rj] 

Degree 2 Degree 3 Degree 4 EE 
Degrees 0-2 0.046 uc SEN = 
Degrees 0-3 0.960 0.794 - E. 
Degrees 0-4 0.884 0.469 0.597 SS 
Degrees 0-5 -0.219 -0.105 -0.035 0051 


The expansion for degrees 0-3 shows the best correlation, but there is no 
consistency from expansion to expansion. Since the least squares solutions do not 
exhibit consistent, stable behavior, it appears that the system Gm = d does not 
provide enough information to uniquely determine the model parameters [or a 


best estimate for them]. This indicates that the system is not even- or 
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overdetermined, but mixed-determined (neither completely overdetermined nor 
completely underdetermined] and requires a more sophisticated method to solve 


for the coefficients. 


Singular Value Decompositon Method 

Theory. Singular value decompositon, or SVD, is one way to solve a 
mixed-determined problem. Its purpose is to partition the system of equations 
into an overdetermined part [that can be solved in the least squares sense] and an 
underdetermined part [that can be solved assuming some a priori information]. 
For the general equation Gm = d, it is like a transformation to the system G'm' 
= d', where m' is composed of an overdetermined part, m? and an 


underdetermined part m" (Menke, 1989): 


О О 
em de > Gum d > 6 o [ms]. d] 


o GU! 


SVD decomposes the data kernel matrix G into three matrices: G = 
UAVÍ. The matrix U is an NxN matrix of orthonormal [orthogonal and of unit 
length] eigenvectors that span the data space S(d). Similarly, the matrix V is an 
MxM matrix of orthonormal vectors that span the model parameter space S(m). 
The matrix A is an NxM diagonal eigenvector matrix with nonnegative diagonal 
elements called singular values, arranged in order of decreasing size. Some of 


the singular values may be zero, making it easy to partition the matrix into a 
submatrix Ap with p nonzero singular values, and several zero matrices: 

Ap O 
70 
where U, and v, are the first p columns of U and V, respectively. 





A -| . This simplifies the data kernel decomposition to: G = UD AVI 
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For the equations Gm = d, the solution is: mest = Ea vg called the 
natural solution (Menke, 1989). If the equation GM = d is to some degree 
underdetermincd, Ap specifies the combinations of model parameters for which 
the equation does provide information; these combinations lie in a subspace of 
the model parameter space 5pm). On the other hand, if GM = d is to some 
degree overdetermined, then Ap specifies the combinations of model parameters 
that the product Gm is capable of resolving; these products span a subspace of 
the data space 5,(4). If none of the singular values arc zero, there arc 
undoubtedly some very close to zero that are affecting the solution variance. 
One way to reduce the solution variance 1s to select a cutoff size for the singular 
values and exclude any singular values smaller than this lie. artificially decide the 


size of p, the number of nonzero singular values]. This is equivalent to throwing 


away some combinations of the model parameters [thus reducing the sizes of Un 
and Yal Hovvever, if the singular values excluded are small, then the solution 
will be close to the natural solution, though the data and model resolution will be 
worse. This is a classic trade off situation between resolution and variance 
(Menke, 1989). 

[t is also possible to dampen the smaller sineular values instead of 
throwing them away [equivalent to the damped least squares method]. The 
drawbacks to this method are that the solution is no longer close to the natural 
solution, the data and model resolution are worse and the damping parameter 
must be determined by trial-and-error. For this study, various methods are used 
to try to determine the optimum number of singular values to keep [p] and all 
singular values with index » p are dropped. 

Desired number of singular values. The first step in determining the 
desired number of singular values is to look at the data kernel spectrums [plots of 


the size of the singular values versus their index] for the mantle component data 
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kernel and the geophysical data kernels (Figs. 3.22-3.25). For the mantle 
components, the data kernel G is only a function of location, so it is the same for 
all four components. For the geophysical data, the data kernels are constructed 
differently, so that all three equations with geoid, gravity and gravity gradient 
data are solving for the same spherical harmonic coefficients. With respect to 


the mantle component data kernel, terms in the geoid, gravity and gravity 


| n. 4 GM (1:1) 
gradient data kernels are multiplied by the additional factors of R, p? an 


- GM (141)(42) 
R 3 


d 


, respectively. 

For comparison, spectrums for the degrees 0-1, 0-2, 0-3, 0-4 and 0-5 
expansions are all plotted, but the emphasis here will be on getting reasonable 
results using the degrees 0-5 expansion. All three geophysical spectrums and the 
geochemical spectrum for this expansion show the singular values graduall y 
decreasing in value, with the last five or so singular values being very close to 
zero. There 1s no obvious cutoff size for the singular values apparent in these 
plots, so other methods must be used to estimate p. | 

For the geophysical control set, it is possible to find the number of 
singular values p needed to most closely approximate the actual coefficients. The 
root mean square error between the actual and estimated geophysical coefficients 


IS given by: 





est 





J 





coefficient rms error — 


where M is the number of coefficients [model parameters]. A plot of coefficient 
rms error versus the number of singular values retained (Fig. 3.26) indicates that 
30, 26 and 14 singular values should be retained, for geoid, gravity and gravity 


gradient, respectively, to most closely approximate the actual coefficients. These 





Íewer singular values it takes for the IMS error to explode [at leas; for these 
Sparse data Sets]. 

Since the coefficients for the gcochemistry data A ана IS NO 
DEMO measure how Closely the estimated coefficients match the actual 
coefficients. What can be done is to (ry to match the Observed data as closely as 
possible, while keeping the Solution variance ál a minimum. As a first step, 
trade-off curves are constructed to bracket the range of p values that balance the 
Size of the model variance and the Spread of the model resolution (Figs. 3.27. 
3.30). The Size of the model variance js based upon the unit covariance matrix 
of the model Parameters, which characterizes the depree of error amplification 
that occurs in the mapping from data to model parameters (Menke, 1989). 
Assuming that the data within the four mantle component vectors and the three 
geophysical Vectors are uncorrelated and have uniform varıance 02 [a reasonable 
‘SSumption for the mantle component vectors based upon the findings in Chapter 
2], the Covariance matrix of the model Parameters is piven by: 

İcov mest) = G cov dG a GG EGET 
Mére GE ¡s the generalized Inverse, which for singular value decomposition IS: 
Gé= VAJ! U/ 
unit covariance Matrix is: 
[cov, mest] = oA cov mest} — G 8G eT _ VA Vi 
ally, the size Of the model varlance is: 


4 M 
size( cov mest/) = É Var, mest E = De [ Var, mest]. = D cov, mesti. 
i=] 


[=] 





e o. U 
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the diagonal elements of the model parameter covariance matrix. With 
increasing values of p, the size of the model variance will increase. 

Since resolution is optimal when the resolution matrices are identity 
matrices, it is possible to quantify the spread of model resolution based on the 


size of the off-diagonal elements of the model resolution matrix R (Menke, 
1989): 


M M 
spread(R) =||R - 15 = 2) 37 İR? - Lil” 
iz] jsl 


where I is the identity matrix and R < V, V], [mest < Rm], With increasing 
values of p, the spread of the model resolution will decrease. 

Trade-off curves of size of model variance versus spread of model 
resolution, as a function of the number of singular values retained, show two 
asymptotes [retaining all 36 singular values gives the largest model variance size] 
(Figs. 3.27-3.30). The ideal range for p, to balance the two measures, is in the 
transition between the asymptotes (Table 3.5). 

Another way to try and pin down the desired number of singular values 
[to most closely approximate the data] is to look at plots of model rms error and 
a variance measure versus the number of singular values retained (Figs. 3.31- 
3.34). Model rms error is given by: 

N 
> (apes qeri 





ə 
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model rms error = 


N 
where 


dP = Gmest = GGEd™s - (U,A, VT (V, AZ! UT) q9*s 2 U,UT q95s 


While V, and UZU, are the identity matrix, V and U,U} are not 


necessarily the identity matrix, since U, and V, do not in general span the 


complete data and model spaces (Menke, 1989). The variance measure used is: 





P 
variance measure = x | AZ, 
te] 


since the solution varıance ıs proportional to A Again , thc goal is to usc thc 
plots of thesc two quantitics to sclect p so that the modcl rms error and the 
solution variance arc balanced (Table 3.5). 

Choosing ranges for p using trade-off curves and the model rms/variance 
Curves 1S a Subjective process. The ranges of values are chosen by eye and there 
is no objective way to select an optirnal value of p from these methods. To make 
the process morc objective, Jacobson and Shaw (1991) suggest applying a 
sequential F-test to SVD problems to find the statistically optimal solution. 
Given a null model with q parameters and a larger general model with b 
parameters [b > q), testing the null hypothesis that the additional [5 - 4] 
parameters in the general model do not improve the fit to the data [compared to 


the null model) requires the use of the F-statistic: 


(RSS, - RSS») (n m h) 


E 
(b - 4) SSA 


where RSS, and RSS) denote the residual sum of squares for the null and general 


models, respectfully, and n is the total number of parameters. F has an F- 
distribution with (b - q,n - b) degrees of freedom. The residual sum of squares 
for a given model is defined as: 


N ^ 
rss =" (gs. are 


tz] 


Values of F can be converted into the probability that the null hypothesis 1s truc, 
ic. that the extra parameters do not result in a better fit. Then the quantity |) 


prob(null hypothesis truc)| is the significance level of the additional parameters. 
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For SVD, the sequential F-test starts by testing the significance of a model 
retaining one singular value against a model retaining no singular values, then 
continues to test models retaining incrementally more singular values against the 
current null model. When a model has reached the 9596 significance level 
[chosen for this application] or higher, it becomes the null model against which 
subsequent models are to be tested, until another model also reaches or surpasses 
95% significance and takes its place. Figures 3.35-3.41 show the F-test results 
for the geophysical and geochemical data sets and Table 3.5 lists the resulting 
optimal p values. In general, it appears that the smoother functions [longer 
wavelength] have higher numbers of significant singular values. 

For determining the value of p, the three different methods agree quite 
well (Table 3.5). The trade-off curves define the largest interval for p, which 15 
constrained further by the model rms/variance curves. For every data set, 
except filtered gravity, the value of p determined by the F-test falls within the 
chosen range of the model rms/variance curves. Even so, the F-test p value for 
filtered gravity does not fall far outside the model rms/variance range |p = 29 
compared to 25] and it does fall within the trade-off range. Since the F-test p 
values are in agreement with the other methods and are by far the most objective 
estimate from the three methods, these values will be used in calculating the 
spherical harmonic coefficients. 

Application. How well the estimated spherical harmonic coefficients of 
the constructed geophysical data sets correlate with the actual GEM-L2 
coefficients is an indicator of how closely the estimated geochemistry coefficients 
may be expected to approximate their true coefficients. Three sets of 
geophysical SVD coefficient solutions are all correlated with the GEM-L2 
coefficients: those that minimize the coefficient rms error and those that 


minimize the model error [selected p values from the F-test] for the filtered and 
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unülrered data sets (Table 3.6). Remember that the data kernel matrices G for 
gravity and gravity gradient are modified so that their spherical harmonic 
coefficients are also estimates of the GEM-L2 coefficients. The correlation 


coefficients r; are calculated as outlined above. Plots of rj versus degree include 


confidence levels based upon a student's t-test. The test statistic for the t-test Is: 


85201 rY 2 


T 


Vı-z di-d 


where n 1s the number of coefficients at that particular degrec [(7 - 2) 2 2/]. T 
has a t-distribution with (n - 2) or 2! degrees of freedom. Given a desired 
significance level and the degrees of freedom, the value of T can be looked up in 


a table. Then the value that r, should have to achieve that significance level can 


be calculated and plotted as confidence levels: 
i 


For the plots of r; versus l, the geophysical coefficients estimated by 


ts 


minimizing the coefficient rms error correlate better than those estimated by 
minimizing the model error and, of those, the unfiltered data set correlates better 
than the filtered data set. All three sets of coefficients correlate well with the 
actual GEM-L2 coefficients at degree 2, except for filtered gravity (Figs. 3.42- 
3.44). In all cases, the geoid coefficient estimates correlate the best. In general, 
gravity and gravity gradient correlate better at even degrees, with the exception 
of the filtered coefficients. For the mantle component coefficients, all this 
implies that the degree 2 coefficients are probably good, but beyond that there is 
no guarantee. Of the four mantle component percentage data sets that are 
expanded, the filtered HIMU data set is unique in that it most closely resembles 
the geoid data set in the variation-distance plots (Figs. 3.13 and 3.16). Thus, 
there is a good possibility that at least the degree 3 coefficients for this data set 


are reasonable as well. 
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Correlation coefficients for the actual GEM-L2 coefficients and tlie 
estimated coefficients cannot be calculated at degrees 0 and | because those 
GEM-L2 coefficients are equal to zero. In contrast, the estimates of these 
coefficients from the constructed geophysical data sets are all positive numbers 
the same order of magnitude as the rest of the estimated coefficients. This 
discrepancy is caused by a sampling bias due to the fact that the oceanic islands 
are all hotspot related and hotspots are associated with geoid highs [Richards et 
al., 1988]; no geoid lows are sampled to balance these highs. It is unclear how 
this bias may affect the estimates of the other coefficients. 

The continuous layer model degree 2 "functions" for the constructed geoid 
data set and the mantle component percentages are reconstructed on a five degree 
grid over the globe from 10<0<170 and -180€«4x180 using the calculated 
coefficients and the appropriate S (Figs. 3.45-3.49). It should be noted 
that the contoured values are not actual geoid anomaly values or component 
percentages, but are deviations from the average [degree 0] geoid anomaly value 
Or component percentage [average constructed geoid = 13.7 m; average filtered 
EMI = 0.27; average EMII = 0.17; average filtered HIMU = 0.31; average DMM 
= 0.25]. For comparison, the actual degree 2 geoid is constructed in the same 
way using the GEM-L2 coefficients [average geoid = 0.0 m] (Fig. 3.50). The 
constructed geoid field agrees well with the actual degree 2 geoid, as already 
indicated by the correlation coefficients. For the mantle components, HIMU 
resembles the actual geoid field with two essentially equatorial highs in 
approximately the same locations; EMI and EMI] also have two highs that 
undulate above and below the equator with a longitudinal shift of ~35° to the east 
with respect to the actual geoid [EMI] has less offset than EMI]; and DMM, with 


its two highs and two lows resembles none of the other degree 2 expansions. 
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Of all the mantle component data sets, filtered HIMU has the best chance 
of getting reasonable values for the degree 3 coefficients. The degrees 2-3 
function for filtered HIMU is reconstructed as before (Fig. 3.51). This can be 


compared to the degrees 2-3 geoid reconstructed from the GEM-L2 coefficients 


(Fig. 3.52). 


SUMMARY 
Viewing the distribution of the OIB reservoir as a continuous layer in the 
mantle and using approximation methods to solve for the spherical harmonic 


coefficients of its expansion reveals the following: 


* The mantle end-member component percentage data have a lot of short 


wavelength energy relative to equally limited geoid, gravity and 


gravity gradient control data sets. 


* With the currently available data, solving for the spherical harmonic 
coefficients is a mixed-determined problem, requiring the use of 


singular value decomposition [S VD] to get viable solutions. 


* The F-test is a simple, objective way to determine the number of 


singular values to retain in SVD for the statistically optimal solution. 


e With the current data coverage, only the degree 2 spherical harmonic 
coefficients can be estimated with a reasonable level of confidence 


using S VD. 
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Continuous layer model degree 2 HIMU closely resembles the degree 2 


geoid. 


Continuous layer model degree 2 EMI and EMII resemble a longitude- 


shifted, undulating degree 2 geoid. 


Continuous layer model degree 2 DMM does not resemble the degree 2 


geoid or the degree 2 expansion of any other mantle component. 
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Table 3.2. Separation of the OIB feature data set into 27 Dupal-type features and 
9 DMM-type features, based upon the percentage of the DMM mantle 
component. 


Dupal-tvpe Features DMM 
Gough 3 2Ü 
Shimada Seamount 6.01 
Tristan de Cunha 6.70 
San Felix/San Ambrosio 8.08 
Walvis Ridge ] 1.80 
Azores 13.74 
Trinidade 14.84 
Cook- Austral Islands 16.22 
Tubuai-Austral Islands 16.93 
New England Seamounts 17.00 
Christmas 1 
St. Helena Бил. 
Kerguelen Plateau 175 
Samoa Islands 18.88 
Comores Archipelago 1937 
Fernando de Noronha 19.80 
Marquesas Archipelago 57: 
Pitcairn 21.57 
Cameroon Line 217: 
Sala Y Gomez 2516 
Society Ridge 1841 
Cape Verde Islands 26.06 
Juan Fernandez Islands 26.63 
Crozet Islands 27.34 
Mascareignes 25.57 
Amsterdam/St. Paul 29.20 
Louisville Seamount Chain 31.34 
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Table 3.2. Continued. 


DMM-type Features 


Balleny 
Cocos 


Marion/Prince Edward 
Galapagos Islands 
Ascension 

Ponape 

Hawaiian Islands 


Iceland 
Nunivak 


%o DMM 
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38.53 


40.88 
41.32 
41.69 
46.14 
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Table 3.4. Change in degree 2 spherical harmonic coefficients for the EMII 
percentage data and the geoid anomaly data as the data sets are expanded to 
progressively higher degrees. 


Expansions 


Geoid 
Degrees 0-2 


Degrees 0-3 
Degrees 0-4 
Degrees 0-5 
EMII 

Degrees 0-2 
Degrees 0-3 
Degrees 0-4 


Degrees 0-5 


AS 


-1.491E-05 
-1.511E-05 
-1.825E-05 


3.946E-05 


-0.065243 
-0.076403 
-0.170987 


2.545222 


A) 


-1.049E-06 
2.88 1E-06 
-7.480E-06 


3.357E-05 


-0.063197 
-0.203182 
-0.221924 


10.864696 


pl 


-4.116E06 
1.55 1E-06 
5050506 


-4.237E-04 


0.269622 
0.466347 
0.785142 


-30.064709 


2 
A? 


1.001 E-05 
6.61 0E-06 
6.044E-06 


-1.436E-04 


0.108911 
0.148947 
0.167714 


-20.805890 


2 
B? 


-6.258E-06 
-8.509E-06 
-1.108E-06 


-3.682E-06 


-0.008334 
0.053770 
0.087842 


-10.644588 
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Table 3.6. Summary of correlation coefficients between the GEM L2 
coefficients and three sets of estimated geophysical coefficients that minimize the 
coefficient rms error and that minimize the model error for filtered and 
unfiltered data sets. 


Data Set p Value Degree 2 Degree 3 Degree 4 Degree 5 


Minimizing coefficient rms error 


geoid 30 0.988 0.909 0.744 0.705 
gravity 26 0.981 0.605 0.734 0.310 
gravity gradient 14 0.926 099 0.457 -0.014 


Minimizing model error (F-test) - unfiltered 


geoid ə 0.988 0.863 0.726 -0.031 
gravity 20 0.891 0.264 0.394 -0.145 
gravity gradient 20 00777 -0.084 0.543 0.085 


Minimizing model error (F-test) - filtered 


geoid 24 0029 0:557 -0.256 0.016 
gravity 29 0.356 0.723 -0.331 -0.333 
gravity gradient 24 0.709 0.334 -0.720 -0.491 
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Fig. 3.1. Variation-distance plot for the EMI mantle component showing the 
range of variation in the component percentage with angular distance between 
the feature locations. To account for the variation requires a minimum sampling 
distance of ~ 14.5” [degree 12 expansion]. 





104 


EMH 


0.45 


0.4 


0.35 


A bsolute Difference 





Angle (deg) 


Fig. 3.2. Variation-distance plot for the EMII mantle component showing the 
range of variation in the component percentage with angular distance between 
the feature locations. To account for the variation requires a minimum sampling 
distance of ~ 39? [degree 4 expansion]. 
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Fig. 3.3. Variation-distance plot for the HIMU mantle component showing the 
range of variation in the component percentage with angular distance between 
the feature locations. To account for the variation requires a minimum sampling 
distance of — 14.5” [degree 12 expansion]. 
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Fig. 3.4. Variation-distance plot for the DMM mantle component showing the 
range of variation in the component percentage with angular distance between 
the feature locations. To account for the variation requires a minimum sampling 
distance of ~ 57° {degree 3 expansion]. 
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Fig. 3.5. Three-dimensional plot of the geographic feature principal component 
data. Axes: X = Z1, Y = Z2, Z = Z3. Symbols: + = features inside the 
DUPAL belt, o = features outside the DUPAL belt. 
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Fig. 3.6. Three-dimensional plot of the geographic feature principal component 
data, with symbols distinguishing percentages of the DMM component. Axes: X 
— Zi, Y < 22, Z — Z3. Symbols: + = <10% DMM, x = 10-20% DMM, 
diamond = 20-30% DMM, square = 30-40% DMM, o = 40-50% DMM, A = 
>50% DMM. Most striking is the large spatial separation in the 30-40% 
Eulegory between: [1] Louisville 2 31.8496, [2] Balleny = 32.17%, [3] Cocos = 
38.53%. 
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Fig. 3.7. Variation-distance plot for the EMI mantle component for the DUPAL 
features only [<32% DMM], showing the range of variation in the component 
percentage with angular distance between the feature locations. Using the 
DUPAL features only shows no reduction in the small-scale variation for the 
EMI component. 
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Fig. 3.8. Variation-distance plot for the EMII mantle component for the 
DUPAL features only [<32% DMM], showing the range of variation in the 
component percentage with angular distance between the feature locations. 
Using the DUPAL features only shows no reduction in the small-scale variation 
for the EMII component. 
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Fig. 3.9. Variation-distance plot for the HIMU mantle component for the 
DUPAL features only [<32% DMM], showing the range of variation in the 
component percentage with angular distance between the feature locations. 
Using the DUPAL features only shows no reduction in the small-scale variation 
for the HIMU component. 
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Fig. 3.10. Variation-distance plot for the DMM mantle component for the 
DUPAL features only [<32% DMM], showing the range of variation in the 
component percentage with angular distance between the feature locations. 
Using the DUPAL features only does show a reduction in the small-scale 
variation for the DMM component, with an increase in minimum sampling 
distance from ~ 57° to 89°. 
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Fig. 3.11. Variation-distance plot for the filtered EMI data set, showing the 
range of variation in the component percentage with angular distance after the 
circular filter is applied. The result is a reduction in the small-scale variation, 
with an increase in minimum sampling distance from — 14.5" to 37? [expansions 


to degrees 12 and 5, respectively]. 
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Fig. 3.12. Variation-distance plot for the filtered EMII data set, showing the 
range of variation in the component percentage with angular distance after the 
circular filter is applied. The result is an increase in the small-scale variation, 
with a decrease in minimum sampling distance from ~ 39° to 27° {expansions to 


degrees 4 and 7, respectively]. 
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Fig. 3.13. Vanation-distance plot for the filtered HIMU data set, showing the 
range of variation in the component percentage with angular distance after the 
circular filter is applied. The result is a dramatic decrease in the small-scale 
variation, with an increase in minimum sampling distance from ~ 14.5? to 83? 
[expansions to degrees 12 and 2, respectively]. 
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Fig. 3.14. Variation-distance plot for the fi 
range of variation in the component percen 
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Fig. 3.16. Variation-distance plot for the constructed geoid data set showing the 
range of variation in the geoid with angular distance between the feature 
locations. To account for the variation requires a minimum sampling distance of 
~ 102° [degree 2 expansion]. 
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Fig. 3.17. Variation-distance plot for the constructed gravity data set showing 
the range of variation in gravity with angular distance between the feature 
locations. To account for the variation requires a minimum sampling distance of 


- 95? [degree 2 expansion). 
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Fig. 3.18. Variation-distance plot for the constructed gravity gradient data set 
showing the range of variation in the gravity gradient with angular distance 
between the feature locations. To account for the variation requires a minimum 
sampling distance of ~ 67° [degree 3 expansion]. 
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Fig. 3.19. Variation-distance plot for the filtered geoid data set, showing the 
range of variation in the geoid with angular distance after the circular filter is 
applied. The filtered geoid data set retains essentially the same angular sampling 
distance [- 93°]. 
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Fig. 3.20. Variation-distance plot for the filtered gravity data set, showing the 
range of variation in gravity with angular distance after the circular filter is 
applied. The filtered gravity data set retains essentially the same angular 
sampling distance [~ 93°]. 
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Fig. 3.21. Variation-distance plot for the filtered gravity gradient data set, 
showing the range of variation in the gravity gradient with angular distance after 
the circular filter is applied. Filtering reduces the small-scale varition in the 
gravity gradient data set, with an increase of angular sampling distance from ~ 
67° to 77° [expansions to degrees 3 and 2, respectively]. 
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Fig. 3.22. Data kernel spectrums for the constructed geoid data kernel G. 
Symbols for the different expansions: - = degrees 0-1, + = degrees 0-2, * = 
degrees 0-3, o = degrees 0-4, x = degrees 0-5. For the degrees 0-5 expansion, 
the singular values approach zéro, but there is no obvious cutoff value. 
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Fig. 3.23. Data kernel spectrums for the constructed gravity data kernel G. 


Symbols for the different expansions: : — degrees 0-1, + = degrees 0-2, * = 
degrees 0-3, o — degrees 0-4, x = degrees 0-5. For the degrees 0-5 expansion, 
the singular values approach zero, but there is no obvious cutoff value. 
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Fig. 3.24. Data kernel spectrums for the constructed gravity gradient data 
kernel G. Symbols for the different expansions: : = degrees 0-1, + = degrees 0- 

= —ideerees 0-3, 0 = degrees 0-4, x = degrees 0-5. For the degrees 0-5 
expansion, the singular values approach zero, but there is no obvious cutoff 
value. 
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Mantle Components Data Kernel Spectrum 
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Fig. 3.25. Data kernel spectrums for the mantle component data kernel G. 
Symbols for the different expansions: : = degrees 0-1, + = degrees 0-2, * = 
degrees 0-3, o = degrees 0-4, x = degrees 0-5. For the degrees 0-5 expansion, 
the singular values approach zero, but there is no obvious cutoff value. 
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Fig. 3.26. Plot of the root mean square error [rms error], as a function of the 
number of singular values retained, between the actual GEM-L2 coefficients and 
those coefficients estimated by the constructed geoid, gravity and gravity 
gradient data sets. Line symbols: — . — . = geoid, - - - - < gravity,'''' — 
gravity gradient, = root mean square of the GEM-L2 coefficients. P values 
minimizing coefficient rms error: geoid = 30, gravity = 26, gravity gradient = 
14. 
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Fig. 3.27. Trade-off curve between model variance and model resolution, as a 
function of the number of singular values retained, for the constructed geoid data 
set. Range for p that balances the two measures is: 15<p<30. Note that trade- 
off curves are determined by the data kernel matrices and so are the same for 


filtered and unfiltered data sets. 


132 


Gravity 
104 , 


103 
102 | 
10! | 
100 | 


10:16 


Size of Model Variance 


102 FR 
103 
104 


10-5 
0 0.1 0.2 0.3 0.4 0.5 0.6 07 0.8 0.9 1 


Spread of Model Resolution 


Fig. 3.28. Trade-off curve betvveen model variance and model resolution, as a 
function of the number of singular values retained, for the constructed gravity 
data set. Range for p that balances the two measures is: 9<p<29. Note that 
trade-off curves are determined by the data kernel matrices and so are the same 
for filtered and unfiltered data sets. 
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Fig. 3.29. Trade-off curve between model variance and model resolution, as a 
function of the number of singular values retained, for the constructed gravity 
gradient data set. Range for p that balances the two measures is: 8<p<26. Note 
that trade-off curves are determined by the data kernel matrices and so are the 
same for filtered and unfiltered data sets. 
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Fig. 3.30. Trade-off curve between model variance and model resolution, as a 
function of the number of singular values retained, for the mantle component 
data set. Range for p that balances the two measures is: 15<p<30. Note that 
trade-off curves are determined by the data kernel matrices and so are the same 
for filtered and unfiltered data sets. 
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Fig. 3.31. Plot of model root mean square error [rms error], as a function of 
the number of singular values retained, between the observed geoid data and the 
geoid data predicted from the calculated coefficients. Balancing the model rms 
error and the model variance gives this range of p: 21<p<25 (filtered and 
unfiltered). Line symbols: = unfiltered model rms error, - - - - = filtered 


model rms error, o—o— = model variance. 
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Fig. 3.32. Plot of model root mean square error [rms error], as a function of 
the number of singular values retained, between the observed gravity data and 
the gravity data predicted from the calculated coefficients. Balancing the model 
rms error and the model variance gives this range of p: 20<p<25 (filtered and 
unfiltered). Line symbols: —— = unfiltered model rms error, - - - - = filtered 
model rms error, o—o— = model variance. 
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Fig. 3.33. Plot of model root mean square error [rms error], as a function of 


the number of singular values retained, between the observed gravity gradient 


data and the gravity gradient data predicted from the calculated coefficients. 
Balancing the model rms error and the model variance gives this range for p: 
20<p<25 (filtered and unfiltered). Line symbols: = unfiltered model rms 
error, - - filtered model rms error, 0—0-— = model variance. 
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Fig. 3.34. Plot of model root mean square error [rms error], as a function of 
the number of singular values retained, between the observed mantle component 
data and the mantle component data predicted from the calculated coefficients. 
Balancing the model rms error and the model variance gives this range for p: 
16€px21 (filtered EMI), 16<p<20 (EMIL, 18xpx23 (filtered HIMU), 19xpx22 
(DMM). Line symbols: = filtered EMI, - - — EMIP s — filtered 
HIMU, - --- < DMM, o—o— = model variance. 
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Fig. 3.35. Plot of F-test significance level as a function of the number of 
singular values retained for the geoid and filtered geoid data sets. Basically, the 
test determines whether additional parameters [singular values] make a 
significant contribution to the model fit of the observed data values. Optimal p 
values [for 95% significance] are: p = 25 [geoid] and p = 24 [filtered geoid]. 
Line symbols: = geoid, - * - * = filtered geoid, - - - - = 95% significance 
level. 
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Fig. 3.36. Plot of F-test significance level as a function of the number of 
singular values retained for the gravity and filtered gravity data sets. Basically, 
the test determines whether additional parameters [singular values] make a 
significant contribution to the model fit of the observed data values. Optimal p 
values [for 95% significance] are: p = 20 [gravity] and p = 29 [filtered gravity]. 
Line symbols: — 5”. = filtered gravity, - - - - = 95% 
significance level. 
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Fig. 3.37. Plot of F-test significance level as a function of the number of 
singular values retained for the gravity gradient and filtered gravity gradient 
data sets. Basically, the test determines whether additional parameters [singular 
values] make a significant contribution to the model fit of the observed data 
values. Optimal p values [for 9596 significance] are: p — 20 [gravity gradient] 
and p = 24 [filtered gravity gradient]. Line symbols: = gravity gradient, 
-:-: — filtered gravity gradient, - - - - = 95% significance level. 
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Fig. 3.38. Plot of F-test significance level as a function of the number of 
singular values retained for the filtered EMI data set. Basically, the test 
determines whether additional parameters [singular values] make a significant 
contribution to the model fit of the observed data values. For filtered EMI, the 
optimal p value [for 95% significance] is: p = 20. Line symbols: o——o - filtered 
EMI, - - - - < 95% significance level. 
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Fig. 3.39. Plot of F-test significance level as a function of the number of 
singular values retained for the EMII data set. Basically, the test determines 
whether additional parameters [singular values] make a significant contribution 
to the model fit of the observed data values. For EMII, the optimal p value [for 
95% significance] is: p = 16. Line symbols: o——o - EMII, - - - - = 95% 
significance level. 
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Fig. 3.40. Plot of F-test significance level as a function of the number of 
singular values retained for the filtered HIMU data set. Basically, the test 
determines whether additional parameters [singular values] make a significant 
contribution to the model fit of the observed data values. For filtered HIMU, the 
optimal p value [for 9596 significance] is: p = 23. Line symbols: o——o z filtered 
HIMU, - - - - = 95% significance level. 
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Fig. 3.41. Plot of F-test significance level as a function of the number of 
singular values retained for the DMM data set. Basically, the test determines 
whether additional parameters [singular values] make a significant contribution 
to the model fit of the observed data values. For DMM, the optimal p value {for 
95% significance] is: p = 22. Line symbols: o—o = DMM, - - - - = 95% 
significance level. 








146 


| Solutions Minimizing Coefficient RMS Error 


Correlation Coefficient 














Degree 


Fig. 3.42. Correlation of geophysics coefficient solutions, that minimize the 

x coefficient rms error with the actual GEM-L2 coefficients. Line symbols: - - - - 
- geoid,....- gravity, - .- . = gravity gradient. Confidence levels are 
determined by a t-test with 2/ degrees of freedom. 
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Fig. 3.43. Correlation of geophysics coefficient solutions, that minimize the 
model rms error for the unfiltered data, with the actual GEM-L2 coefficients. 
Line symbols: - - - - = geoid,.... = gravity, - .- . = gravity gradient. 
Confidence levels are determined by a t-test with 2/ degrees of freedom. 
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Fig. 3.44. Correlation of geophysics coefficient solutions, that minimize the 
model rms error for the filtered data, with the actual GEM-L2 coefficients. 
Line symbols: - - - - = geoid,.... = gravity, - .- . — gravity gradient. 
Confidence levels are determined by a t-test with 2/ degrees of freedom. 
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CHAPTER 4 
SPHERICAL HARMONIC REPRESENTATION OF ISOTOPIC 
SIGNATURES: THE DELTA-FUNCTION MODEL 
INTRODUCTION 

As mentioned in Chapter 3, the delta-function model represents the OIB 
reservoir as a series of point sources, each feeding a separate plume. This may 
seem unphysical, but could be a good approximation of actual conditions if the 
source boundary layer is not continous, but patchy, as indicated in some seismic 
studies of D" (Lay et al., 1990). 

Representing the geographic features as delta-functions [scaled by the 
corresponding geoid anomaly or mantle component percentage] has two 
advantages, mathematically, over the approximation methods used in Chapter 3. 
First, the spherical harmonic coefficients can be found easily with the 
simplification from integration over the globe to summation over the feature 
locations allowed by the delta-functions. Second, representing the OIB reservoir 
as a known function removes the problem of aliasing; the values of the spherical 
harmonic coefficients are not dependent upon the truncation point of the 
expansion [they are dependent upon the number and location of the geographic 
features]. For delta-functions, which have energy at all degrees, the expansions 
can be carried out to infinity, but for this study, will only be carried out to 


degree 5, for comparison with the continuous layer model. 


THEORY 


As before, any function f(8,@) can be expanded in spherical harmonics: 


5 
END Al m P/'(cos0) [ A7'cosmqQ 4 B7'sinmg] 


120) m=0 





TOB 


Due to the orthogonality of the spherical harmonics, the equations for the 


coefficients are: 


AT - m p) q EA e P7 (cos0) cosmp d(cos0) 
Ar(I+m)! 


Bj -| el f( 0.) A / y" 3 Р/"(со50) sinn d(cos0) 


For the delta-function model, the function being expanded is a series of delta- 


functions; 


R0,9) = k;0-0,,9-9,) 


where X; is one of the four mantle component percentages [or the value of the 


geoid anomaly] and X6-0;,0-0;) indicates a delta-function at the particular 


location (8;,9;). Mathematically, the delta-function is a “spike” of infinite height, 


infinitesimal width and unit area: 


| 2 ö(0-0),9-pi) dö < 1 


The key property of the delta-function is that the integral of a function 2(6,q@) 


times a delta-function is just the value of 2 at the delta-function location: 


| аф J £(0,9) 0(0-0;,9-«;) dO — 9(0/,;) 


This simplifies the coefficient equations from integration over the globe to 


summation over the geographic feature locations: 


A = y 7: oU Pj (cos0;) cosmt; 


i=l 





KOG 


N 205160 
5 ES k, C P7"(cos0;) sinmo, 


The coefficient equations for the constructed data sets of geoid, gravity 
and gravity gradient anomalies at the feature locations have additional factors. 


As an example, for gravity the equations are: 


(21-1) -m)! 
Ay = A Pr(cos0;) cosmy; 
TUUS je aiit a) Ee, 
(2-1)- m)! 
Bj LÀ. PI (cos0;) sinm 
' "sihi N dm Lr) (cose) É 
R? 
with the additional factor of GM(l+1). Geoid and gravity gradient additional 
R? 


] һжж— 
factors are p and GM(l+1X1+2), respectively. 


APPLICATION 

As before, the constructed geophysics data sets are used as a control to 
gauge the level of accuracy expected from the mantle component data sets. 
Correlating the coefficients from these data sets with the GEM-L2 cocificients 
(Fig. 4.1) yields good agreement for all three at degree 2. Whereas the 
continuous layer model showed a fairly consistent pattern of decreasing 
correlation from the geoid coefficient estimates to the gravity and gravity 
gradient estimates (Fig. 3.43), the delta-function model shows equal correlation 
at degree 2 and a switch to increasing correlation from the geoid estimates to the 


gravity and gravity gradient estimates at degree 4. Overall, it appears that the 





delta-function model is less accurate at reproducing the coefficients for long 
wavelength data sets [geoid] and more accurate at reproducing the coefficients 
for the short wavelength data sets [gravity gradient) than the continuous layer 
model. Both models are consistent, though, in showing strong correlation for all 
three data sets at degree 2, implying that the mantle component degree 2 
coefficients are also viable. In addition, the mantle component data sets have 
even more high degree [short wavelength] energy than the gravity gradient data 
set, so their coefficients are probably reasonably accurate out to degree 4. 

Since each of the different geophysics data sets approximate the GEM-L2 
coefficients equally well at degree 2, it appears that there is some additional 
controlling factor affecting the estimates of the degree 2 coefficients, aside from 
the data values themselves. The location of the features, and thus the delta- 
functions, is the most likely candidate. A plot of the constructed degree 2 
“function” for the delta-function model geoid (Fig. 4.2) shows the obvious 
relationship between the two main clusters of oceanic islands and the two highs in 
the geoid. Since the continuous layer model geophysics coefficients all agreed 
well with the degree 2 GEM-L2 coefficients, it appears that the location effect 
merely enhances an already existing correlation and is not solely responsible for 
the correlation. Presumably the same is true of any degree 2 correlation of 
delta-function model geochemistry coefficients with the GEM-L2 coefficients. 

Degree 2 "functions" for the mantle component percentages are 
reconstructed, as before, for comparison with those of the continuous layer 
model (Figs. 4.3-4.6). The contoured values of the delta-fuction geoid (Fig. 4.2) 
and the mantle component functions are large enough to be the actual geoid and 
component percentages, instead of deviations from the average values, as for the 


continuous layer model. This is due to the arbitrary scaling that comes into play 


ıı 


when using delta-functions. A delta-function has unit area, so the average valuc 
of a delta-function over the sphere 1s: 


= A 
(Ap sin0) ÃO 4r 


where (Aq sin0) A0 is a sectional area on the sphere (Fig. 4.7), which for the 
whole sphere is 47. If there is only one delta-function involved in the 
reconstruction, the contoured values will be off by a factor of 1/(4m). Since 
there are 36 features, there are 36 delta-functions involved in the reconstruction, 
so the contoured values are off by a factor of 36/(47) = 2.86 or ~3. 

Qualitatively, the four reconstructed mantle component degree 2 functions 
show good agreement with each other. All four have two highs: one over central 
Africa and the other over the central Pacific. Slight differences include the 
width of the highs [from narrowest to widest width: HIMU, EMI, EMI! and 
DMM] and the amount of displacement [from 0° to 15°] of the highs above and 
below the equator [from least to most displacement: HIMU, EMII, EMI and 
DMM]. With respect to the GEM-L2 degree 2 geoid (Fig. 3.50), all of the 
mantle component highs are shifted longitudinally to the east by varying amounts 
[HIMU ~30°, EMI -30°, EMII -35° and DMM -40?]. 

Degrees 2-3 functions for the four components (Figs. 4.8-4.1 1) are 
constructed for comparison with the geoid (Fig. 5.32) and the HIMU continuous 


layer model reconstruction (Fig. 5.31). 


SUMMARY 
Viewing the distribution of the OIB reservoir as a series of point sources 


that can be represented as delta-functions yields the following results: 





With respect to the behavior of geophysics control data sets, at least the 
degree 2 spherical harmonic coefficients for the mantle components 


can be estimated with confidence, if not the degrees 3 and 4 as well. 


The location of the features, and thus the delta-functions, biases the 
calculated degree 2 coefficients due to the correlation between the 


oceanic island locations and the degree 2 geoid. 


Scaling of delta-function models reconstructed over the globe is 
dependent upon the number of delta-functions used in the 


approximation [A] and varies as N/(4z). 


Degree 2 HIMU, EMI, EMIT and DMM all show a degree 2 geoid 
pattern phase-shifted 30°-40° to the east, with varving widths of the 


highs and displacements from the equator. 





Delta-Function Model Geophvsics Correlation with GEM-L2 
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Fig. 4.1. Correlation of the delta-function model geophysics coefficient solutions 
with the actual GEM-L2 coefficients. Line symbols: -- - - = geoid..... = 
gravity, - .- . = gravity gradient. Confidence leveis are determined by a t-test 
with 2/ degrees of freedom. 
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is longitude. 


Fig. 4.7. Geometry of a sectional area on a sphere, where 0 is colatitude and Q 
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CHAPTER 5 


RESULTS AND DISCUSSION 


INTRODUCTION 

Geophysical control data sets are used to judge the dependability of 
spherical harmonic coefficient solutions for the mantle end-member components 
from the continuous layer and the delta-function models. A careful comparison 
of the two models can further enhance or reduce the significance assigned to the 
various solutions. [n this chapter, the two models are compared in terms of their 
amplitude spectra, how well they correlate with the geoid, how they are affected 
by nonuniform feature distribution and how well they correlate with the 
Clayton-Comer seismic tomography model. The significance of the correlations 
with the geoid and the seismic tomography model is discussed, along with 


suggestions for further research. 


AMPLITUDE SPECTRA 
Spectral amplitude plots show the relative power at cach degree for the 
different mantle component expansions. Following Richards and Hager (1988), 


the root mean square harmonic coefficient amplitude at each degree is given by: 


| «(pr 


“Ar. "EL — RH S 


2 7 
where V/ is the variance at each degree for a given set of harmonic coefficients. 





Richards and Hager (1988) include the factor of 1/(2/ 4 1) because random noise 
on a sphere will have a flat spectrum with this normalization. On plots of 
SI versus l, low-degree or long-wavelength effects will show up as a negative 


slope. 





OO 


Amplitude spectra of the calculated geoid coefficients from the two 
models agree well with the negative [long-wavelength] slope of the actuat geoid 
coefficients (Figs. 5.1 and 5.2). For the mantle component expansions, 
amplitude spectra reveal no such clear cut negative slope pattern to indicate 
dominant long-wavelength effects (Figs. 5.3 and 5.4). Instead, the spectra appear 
"white", with energy at all degrees, and no decrease in the energy with 
increasing degree. In addition, HIMU is the only mantle component that shows 
any consistency in behavior between the two models. Thus, in general, the 


expansion of the mantle components is model dependent. 


CORRELATION WITH THE GEOID 

Plotting the mantle component percentages point by point against the full 
geoid value at the geographic feature locations 1s not a valid way to compare the 
mantle component signatures with the geoid. When correlating them by degree 
using spherical harmonic coefficients, it is apparent that the mantle components 
may correlate with the geoid at some degrees [wavelengths] and not others. In a 
pointwise comparison, the different patterns at the different degrees are 
obscured as they are added together to produce the whole, making an accurate 
comparison impossible. Pointwise plots done with the current data show no 
correlation between the mantle components and the geoid (Figs. 5.5-5.8). 

[n contrast, correiating the geoid coefficients and the mantle component 
coefficients by degree reveals a good corrrelation [90% significance level and 
higher] at degree 2 for the DUPAL components [EMI, EMH and HIMU] for both 
models (Figs. 5.9 and 5.10). Note that positive correlations indicate high 
concentrations of mantle components correlating with geoid highs and vice versa. 
HIMU has the best correlation for both models, showing better than 95% 


significance at degree 2 and 90% sienificance at degree 3. The remaining mantle 
E g E - 


İn” 


components show a consistent decreasing correlation from EMII to EMI to 


DMM for both models. 


IMPLICATIONS OF NONUNIFORM FEATURE DISTRIBUTION 

Oceanic island distribution is not uniform about the globe. As indicated in 
Chapter 4, the two main clusters of oceanic islands correspond to the two highs 
of the degree 2 geoid. It can be argued, then, that any correlation between the 
degree 2 mantle component expansions and the degree 2 geoid is due solely to the 
nonuniform distribution of the oceanic islands and not to any pattern in the 
geochemistry values. Fo test this, the percentages of the HIMU mantle 
component at the 36 geographic features, filtered [continuous layer model| and 
unfiltered [delta-function model], are randomly assigned to different feature 
locations five times. HIMU percentages are used since the degree 2 HIMU, for 
both models, correlates best with the degree 2 geoid. The five randomly 
generated data sets for each model are then used to compute new coefficients that 
can be compared to the degree 2 geoid. For the continuous layer model, the 
number of singular values retained for the new data sets is determined by the F- 
test at the 95% significance level. The random number generator used for this 
test is nonlinear, but repeatable, since it starts with a given seed that is updated 
for successive calls in a predictable manner. This means that for a given 
randomization, the filtered and unfiltered HIMU percentages are being 
randomized in the same way, so the results of the two models can be compared. 
Five iterations is not enough to quantify the effect of the feature distribution on 
the degree 2 correlation for the two models, but it is enough to indicate if it has 
any control at all. 

Concentrating on the degree 2 coefficients, three of the randomizations 


that result in strong correlations with the geoid for delta-function model [well 
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above the 90% confidence level] result in negligible correlations with the gcoid 
for the continuous layer model (Table 5.1). Reconstructed degree 2 functions of 
the randomized data sets show graphically how little the delta-function model 
changes, with respect to the continuous layer model, when the geochemical 
signatures of the features are mixed up (Figs. 5.11-5.20). For the delta-function 
model, this indicates that the values of coefficients are not so much dependent 
upon the scaling factors multiplying the delta-functions as the location of the 
delta-functions themselves. This location effect makes it difficult to trust strong 
correlations of the delta-function model with the geoid unless there 1s additional 


confirmation by the continuous layer model. 


CORRELATION WITH SEISMIC TOMOGRAPHY 

Correlating the mantle component expansions with the geoid gives an 
estimate of the general OIB source region [ie. lower mantle versus upper 
mantle], but is incapable of resolving a more precise depth range for the source 
since the geoid 1s affected by mass anomalies at all depths in the Earth. A way to 
select a probable depth range for the OIB source[s] is to compare the mantle 
component expansions to seismic tomography models. Seismic tomography 
models map the global distribution of lateral velocity variations in the mantle at 
different depths based upon the inversion of travel time anomaly data from 
seismic waves that travel through the Earth's interior (Hager and Clayton, 1989). 

In this study, the mantle component expansions are correlated with the 
Clayton- Comer seismic tomography model, discussed in Hager and Clayton 
(1989). The Clayton-Comer model inverts for slowness [inverse of velocity] 
anomalies, in a given shell, that are converted to velocity anomalies by 
multiplying by the average shell velocity. There are 29 shells in the model, each 


100 km thick, spanning the entire mantle from the core-mantle boundary [CMB], 
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at a depth of 2900 km, to the surface. Shells 23-29 [covering the uppper mantle] 
are not used in this analysis since coverage in the top 700 km of the mantle is 
poor because of the near vertical seismic ray paths in this region. The spherical 
harmonic coefficients of the remaining 22 shells [covering the lower mantle] are 
averaged together, to dampen model noise, to produce 5 layers: 2900-2500 km 
[layer 1], 2500-2100 km [layer 2], 2100-1700 km [layer 3]. 1700-1200 km [layer 
4] and 1200-700 km [layer 5]. 

The geoid is correlated with the Clayton-Comer tomography model first 
(Fig. 5.21) to serve as a guide for interpreting the correlation of the tomography 
model with the mantle component expansions. Note that a negative correlation 
indicates geoid highs correlating with low velocity regions [and vice versa] and a 
positive correlation indicates geoid highs correlating with high velocity regions 
[and vice versa]. In layers 1-3, the strong negative correlations at degrees 2 and 
3 confirm that long wavelength geoid highs are due to low density [warmer and 
thus slower velocity] mantle upwellings. This long wavelength upwelling 
signature is also present in the upper lower mantle, as shown by the strong 
negative correlations at degrees 2, 3 and 4 for layer 4 and at degree 2 for layer 
5. Of interest is the strong positive correlations for layers 4 and 5, at degree 5 
and degrees 4 and 5, respectively. Bowin (19913) indicates the correspondence 
of the degrees 4-10 geoid highs with plate convergence zones. He believes that 
the mass anomalies responsible for the highs lie in the lower mantle, beneath 
plate convergence zones, below the teleseismically downgoing subducted slabs. 
The positive correlations in layers 4 and 5 support this theory and imply that 
subducted slabs extend below the 670 km discontinuity. 

Correlation of the mantle component expansions with the Clayton-Comer 
tomography layers for the two models yields interesting results (Figs. 5.22- 


5.29). Due to the limitations of both models [ie. the uncertainties in the 


——— — ... 
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coefficient estimates for the continuous layer model and the location dependence 
in the delta-function model], it is more likely that a significant correlation is 
accurate if it is present in both models. With this in mind, the interpretation of 
the correlation results will be based upon common correlations of 90% 
significance [or very close to it] or higher (Table 5.2). 

The common degree 2 correlations with layers 3-5 for all the mantle 
components are indicative of large scale upwelling, as for the geoid. Good 
degree 3 correlations with layer 1 points to a deep source for all four 
components, like the geoid which shows a much stronger correlation at degrec 3 
with layer ] than 1t does at degree 2. This correlation is not unexpected for the 
DUPAL components, whose correlation with the degree 2 geoid also suggest a 
deep origin, but it is surprising for the DMM component. There are two 
possible solutions for the dilemma posed by the supposedly upper mantle DMM 
component correlating with deep mantle tomography. First, it is possible that 
the DMM component expansion does correlate better with upper mantle 
tomography, which is, unfortunately, not available for the Clayton-Comer 
model. Second, it is possible that the DMM component is representative of both 
the upper and lower mantle composition. Hart (1991) shows that all the hotspots 
that have elongated isotopic arrays indicate mixing between one of the DUPAL 
components and something that is nota MORB composition. Since 3/4 of a 
plume’s ascent is spent in the lower mantle, the composition of the DMM 
component may be largely controlled by lower mantle entrainment (Hart, 1991). 

Another interesting correlation common to both models is the positive 
correlation at degree 5 for EMIT in layer 5. With respect to Bowen's model 
(19912) this indicates a correlation between the EMII component and subducted 
slabs. This finding agrees with the geochemical evidence suggesting the EMII 


component is derived from recycling of subducted sediments (Hart, 1988). 
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DISCUSSION 

As indicated in Chapter 3, the average value of the geoid anomaly at the 
36 feature locations is 13.7 m, not zero as it should be if the features were 
located randomly with respect to the geoid. Fhis 1s a simple indication that the 
feature locations [hotspots] correlate with geoid highs. Naturally, then, the bulk 
chemical signatures unique to oceanic island basalts should also correlate with 
geoid highs. What is significant is that the expansions of all three DUPAL 
mantle end- member components [EMI, EMH and HIMUJ], that comprise 3/4 of 
the bulk chemical signature, individually correlate with geoid highs. More 
importantly, the DUPAL components correlate with the degree 2 geoid highs, 
indicating a deep origin for the components since the degrees 2-3 geoid ficld is 
inferred to result from topography at the core-mantle boundary (Bowen, 19914). 

It can be argued that the correlation of the DUPAL components with the 
degree 2 geoid is not an indication of geochemical patterns within the earth, but a 
direct result of the nonuniform distribution of the oceanic islands, whose two 
largest population densities correspond to the degree 2 geoid highs. 
Randomization tests indicate, however, that while this nonuniform distribution 
does play a role in solutions for the delta-function model, it 1s not the controlling 
factor for continuous layer model solutions. Though the continuous layer model 
solutions are hindered by the limited number and coverage of the oceanic islands 
and the delta-function model solutions are biased by the oceanic island locations, 
continual comparisons of the two models can be used to judge the accuracy of the 
solutions [in addition to judging accuracy using geophysical control sets]. 
Essentially, where both models agree, the solutions are more likely to be 


accurate. 
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The total geoid field is due to the contribution of different mass anornalies 
at different depths throughout the Earth, so it can be difficult to directly 
ascertain a source depth by comparing geochemical quantities with the geoid. 
Seismic tomography models allow the correlation of geochemical quantities with 
seismic velocity anomalies at different depths and serve as an independent check 
on the general source locations indicated by correlation with geoid anomalies. 
Correlating the mantle end-member components from both models with the 
Clayton-Comer seismic tomography model suggests a source depth range of 
2500-2900 km [just above the core-mantle boundary] for the DUPAL 
components, dus to the strong negative degree 3 correlations at this depth. In 
addition, a strong positive degree 5 correlation in the depth range of 700-1200 
km is an indication that the EMJI component is related to subduction, as 
previously suggested using geochemical evidence (Hart, 198%). Sirnilarly, the 
geoid shows a strong positive correlation with the Clayton-Comer model at 
degrees 4 and 5 in the depth range 709-1200 krn and at degree 5 in the depth 
ranges of 1200-1700 km. These subduction related patterns in the upper lower 
mantle indicate that subducted slabs extend beyond the 670 km seisrnic 
discontinuity and thus are supporting evidence for whole mantle convection 

Further comparisons need to be rnade between the mantle component 
expansions and other seisrnic tomography models. ft is especially important to 
compare the mantie components to a high resolution upper mantle tomography 
model, since the arnplitude spectra for the cornponents indicate power at high 
degrees which will becorne dorninant at shallow depths in the mantle. Such a 
comparision could clarify the nature of the DMM component, which correlates 
well with the degree 3 deep mantle layer of the Clayton-Comer model, and could 


further explore the relationship between the EMIT component and subduction. 
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SUMMARY 
A comparison of the two models used to expand the mantle components in 


spherical harmonics yields the following results: 


* Mantle end-member component amplitude spectra, for the continuous 
layer model and the delta-function model, show power at all degrees, 


with no one degree dominating. 


. The DUPAL components [EMI, EMII and HIMU] for both models 


correlate well with the geoid at degree 2, indicating a deep origin. 


e  Delta-function model solutions are, to some extent, controlled by the 
nonuniform feature distribution, while the continuous layer model 


solutions are not. 


* The DUPAL and DMM components, for both models, correlate well 
[negatively] at degree 3 with the velocity anomalies of the Clayton- 
Comer seismic tomography model in the 2500-2900 km depth range 


[immediately above the core-mantle boundary]. 


* The EMII component, for both models, correlates well [positively] at 
degree 5 with the velocity anomalies of the Clayton-Comer seismic 
tomography model in the 700-1200 km depth range, indicating a 


subduction related origin. 


° Subduction related positive correlations for the geoid and the EMII 


component with the Clayton-Comer model in the upper lower mantle 
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[700-1700 km] indicate that subducted slabs extend below the 670 km 


seismic discontinuity, Supporting a whole-mantle convection model. 
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able 5.1. Summary of correlation coefficients between the GEM-L2 
oefficients and coetticients calculated from five randomly generated data sets 

r the continuous laver model [filtered HIMU] and the delta-function model 
IIMU], along with the actual correlations of the filtered HIMU and HIMU data 














(S. 

lata Set Degree 2 Degree 3 Degree 4 Degree 5 
ontinuous Layer Model 

Itered HIMU 0557 0.502 -0.112 -0.358 
ındom | 0.753 0.446 -0.639 -0.157 
indom 2 0.560 0.196 0.467 -0.210 
indom 3 -0.129 0.432 -0.303 -0.069 
indom 4 0223 0.448 0.386 -0.130 
indom 5 -0.166 0.718 -0.285 -0.230 
Jelta-Function Model 

IMU 0.850 0.491 0.063 -0.505 
andom 1 0.726 0.332 0.036 -0.416 
andom 2 0.622 0.404 0.107 -0.361 
andom 3 0.873 0.477 0.029 -0.320 
andom 4 0.893 0.407 0.415 -0.385 
indom 5 0.761 0.383 0.107 -0.286 


_—_ 0807 — ———————————— 
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Table 5.2. Summary of correlations of 90% significance [or very close to it] or 
higher for the continuous layer model and the delta-function model when 
correlated with five averaged layers in the Clayton-Comer tomography model. 
[A "+" or "-" next to the component name indicates a positive or negative 
correlation, respectivelv.] 


Degree 2 Degree 3 Degree 4' "Pesrce 5 


-EMI 
Laver 5 -EMII +EMII 
-HIMU! 
-EMI 
Laver 4 -EMII! 
-DMMİ 
-EMI 
Laver 3 -EMII 
Laver 2 -DMM? -EMIIT 
-EMI! -EMI 
Layer 1 -EMII! -EMIT 
-HIMU 
-DMM -DMM 


The continuous layer model correlation is slightly less than 90% significant. 
2The delta-function model correlation is slightly less than 90% significant. 
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x10:6 Continuous Layer Model Geoid 


RMS Amplitude 
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Fig. 5.1. Amplitude spectra for the continuous layer model coefficients of the 
constructed geoid data set, as compared to the actual geoid. Line symbols: - - - - 
= GEM-L2 geoid. 





= constructed geoid, 
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Delta-Function Model Geoid 
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Fig. 5.2. Amplitude spectra for the delta-function model coefficients of the 
constructed geoid data set, as compared to the actual geoid. Line symbols: - - - - 


= constructed geold, 





= GEM-L2 geoid. 
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Conunuous Layer Model Mantle Components 
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Fig. 5.3. Amplitude spectra for the continuous layer model coefficients of the 
mantle component data sets. Line symbols: = filtered EMI, -- - - = 
EMIL - < fihered HIMU, - ' -: < DMM. 
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Delta-Function Model Mantle Components 
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Fig. 5.4. Amplitude spectra for the delta-function model coefficients of the 
mantle component data sets. Line symbols: = EMI, ----=EMII,----= 


HIMU, - --- = DMM. 
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Fig. 5.5. Pointwise comparison, at each geographic feature, of the full geoid 


anomaly [in meters] with the EMI component percentage. This plot gives the 
impression that there is no correlation. 
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Fig. 5.6. Pointwise comparison, at each geographic feature, of the full geoid 
anomaly [in meters] with the EMII component percentage. This plot gives the 
impression that there is no correlation. 
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Geold 





Fig. 5.7. Pointwise comparison, at each geographic feature, of the full geoid 


anomaly [in meters] with the HIMU component percentage. This plot gives the 
impression that there is no correlation. 








Fig. 5.8. Pointwise comparison, at each geographic feature, of the full geoid 
anomaly [in meters] with the DMM component percentage. This plot gives the 
impression that there 1s no correlation. 





Correlation Coefficient 


Fig. 5.9. Correlation of the continuous layer model mantle component 
coefficient solutions with the GEM-L2 geoid coefficients. Line symbols: 
filtered EMI, - - - - = EMII, :: +: = filtered HIMU, ---- = DMM. Confidence 
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Continuous Layer Model Mantle Components Correlation with GEM-L2 
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levels are determined by a t-test with 2/ degrees of freedom. 
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Delta-Function Model Mantle Components Correlation with GEM-L2 


Correlation Coefficient 





Fig. 5.10. Correlation of the delta-function model mantle component coefficient 
solutions with the GEM-L2 geoid coefficients. Line symbols: — EMI, --- 
- EMI,;::: 2 HIMU,-:--  DMM. Confidence levels are determined by a 
t-test with 2/ degrees of freedom. 
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Correlation of Geoid vvith Clayton-Comer Model 
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Fig. 5.21. Correlation of the GEM-L2 geoid coefficients with the five layers of 
the Clayton-Comer seismic tomography model. Line symbols: Slayer l 
[2500-2900 km], - - - - = layer 2 [2100-2500 km], : +++ = layer 3 [1700-2100 
km], - - - - = layer 4 [1200-1700 Km jo o jal co 1700-1200 km]. 
Confidence levels are determined by a t-test with 2/ degrees of freedom. 
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Fig. 5.22. Correlation of the continuous layer model filtered EMI coefficients 
with the five layers of the Clayton-Comer seismic tomography model. Line 
symbols: = layer 1 [2500-2900 km], - - - - = layer 2 [2100-2500 km], + -> 
- = layer 3 [1700-2100 kmj, - ' - : < layer 4 (1200-1700 km], o—o = layer 5 
[700-1200 km]. Confidence levels are determined by a t-test with 2/ degrees of 
freedom. 
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Continuous Layer Model EMII Correlated with Clayton-Comer Model 
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Fig. 5.23. Correlation of the continuous layer model EMII coefficients with the 
five layers of the Clayton-Comer seismic tomography model. Line symbols: 

= layer 1 [2500-2900 km], - - - - = layer 2 [2100-2500 km], - - - - = layer 
3 [1700-2100 km], - - - : = layer 4 [1200-1700 km], o—o = layer 5 [700-1200 
km]. Confidence levels are determined by a t-test with 2/ degrees of freedom. 
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Continuous Layer Model HIMU Correlated vvith Clayton-Comer Model 
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Fig. 5.24. Correlation of the continuous layer model filtered HIMU coefficients 
with the five layers of the Clayton-Comer seismic tomography model. Line 
symbols: = aver 1 (25002900 km], == layer 2 [2100-2500 km], - : : 
aer 3 [1700-2100 kml, - * -- = layer 4 [1200-1700 kmi, 6—— o z layer > 
(700-1200 kml. Confidence levels are determined by a r-test with 2/ degrees of 
freedom. 
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Continuous Layer Model DMM Correlated vvith Clayton-Comer Model 
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Fig. 5.25. Correlation of the continuous layer model DMM coefficients with the 
five layers of the Clayton-Comer seismic tomography model. Line symbols: 

= layer 1 12500-2900 kmi, - - - - — layer 2 [2100-2500 km], ---- = layer 
3 [1700-2100 km], - : -  = layer 4 [1200-1700 km], o—o layer 5 1700-1200 
km]. Confidence levels are determined by a t-test with 2/ degrees of freedom. 
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Delta-Function Model EMI Correlated with Clayton-Comer Model 
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Fig. 5.26. Correlation of the delta-function model EMI coefficients with the five 
layers of the Clayton-Comer seismic tomography model. Line symbols: = 
layer 1 [2500-2900 km], - - - - < layer 2 [2100-2500 km], ': : : = layer 3 [1700- 
2100 km], - : - : — layer 4 (1200-1700 kmi, o-—o = layer 5 1700-1200 kml. 
Confidence levels are determined by a t-test with 2/ degrees of freedom. 
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Delta-Function Model EMII Correlated with Clayton-Comer Model 
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Fig. 5.27. Correlation of the delta-function model EMII coefficients with the 
five layers of the Clayton-Comer seismic tomography model. Line symbols: 

= layer 1 [2500-2900 km], - - - - = layer 2 [2100-2500 km], : = layer 
3 [1700-2100 km], - : - : = layer 4 [1200-1700 km], o——o = layer 5 [700-1200 
km]. Confidence levels are determined by a t-test with 2/ degrees of freedom. 
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Da Delta-Function Model HIMU Correlated with Clayton-Comer Model 
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Fig. 5.28. Correlation of the delta-function model HIMU coefficients with the 
five layers of the Clayton- Comer seismic tomography model. Line symbols: 
——— z laver 1 [2500-2900 km], - - - - = layer 2 [2100-2500 km], : : :: = layer 
3 [1700-2100 km], - : - : = layer 4 [1200-1700 km], o——o = layer 5 [700-1200 
km]. Confidence levels are determined by a r-test with 2/ degrees of freedom. 
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Delta-Function Model DMM Correlated vvith Clayton-Comer Model 
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Fig. 5.29. Correlation of the delta-function model DMM coefficients with the 
five layers of the Clayton-Comer seismic tomography model. Line symbols: 

= layer 1 [2500-2900 km], - - - - = layer 2 [2100-2500 km), >+ >+ = layer 
3 (1700-2100 km], - : - - 2 layer 4 (1200-1700 km], o——o = :ayer 5 1700-1200 
km]. Confidence levels are determined by a t-test with 2/ degrees of freedom. 
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